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ABSTRACT 

The  stability  of  plane  Poiseuille  flow  was  studied 
using  theory  developed  by  Harrison.   A  similarity  trans- 
formation was  introduced  which  reduces  computation  time 
and  provides  better  insight  into  the  basic  relations. 
The  stability  of  the  flow  was  examined  from  a  Lagrangian 
viewpoint.   Instability  was  found  to  be  progressive  in 
nature  and  three  distinct  levels  were  identified,  namely 
incipient,  critical,  and  fully  developed  instability. 

Results  show  that  the  critical  Reynolds  number  can  be 
lowered  indefinitely  if  certain  types  of  perturbations 
occur.   Specifically  these  involve  relatively  abrupt 
changes  in  amplitude.   This  provides  a  possible  expla- 
nation for  the  disagreement  between  earlier  theory  and 
experiment. 
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LIST  OF  SYMBOLS 

All  quantities  are  expressed  in  dimensionless  form  by 
the  use  of  a  natural  system  of  consistent  units  in  which 
channel  semi-height  is  the  unit  of  length,  the  volumetric 
mean  velocity  of  the  fluid  is  the  unit  of  velocity,  and 
the  density  of  the  fluid  is  the  unit  of  density.   Then 
all  other  consistent  derived  units  are  fixed  accordingly. 

A*  wave  number  amplitude  defined  in  Eq.  2.4 

e  2. 71828. . .base  of  the  natural  logarithms. 

G  stability  parameter  defined  in  Eq.  1-30. 

G,H  complex  stream  functions. 

G*,H*  transformed  stream  functions. 

I (y)        complex  auxiliary  function  defined  in  Eq.  2-14 

1/2 
i  (-1)     the  imaginary  unit. 

r,  J",  E"     unit  vectors  along  x,  y,  and  z  axes, 
respectively. 

J(y)        complex  auxiliary  function  defined  in  Eq.  2-15 

J*(y)       transformed  auxiliary  function  defined  in 
Eq.  2-27. 

Re  Reynolds  number  based  on  volumetric  mean 

velocity  and  channel  semi-height. 

Re*         transformed  Reynolds  number  defined  in 
Eq.  2-18. 

T(y)        complex  auxiliary  function  defined  in 
Eq.  2-13. 

t  time. 

U  flow  velocity. 


W  complex  vector  potential  of  perturbation 

flow  defined  in  Eq.  3-2. 

x,y,z       coordinates  in  direction  of  mean  flow, 

normal  to  walls  and  transverse  to  the  mean 
flow,  respectively. 

x',y,z      coordinates  in  moving  reference  frame. 

a  complex  wave  number  of  the  perturbation  in 

x  direction. 

a*  transformed  wave  number  defined  in  Eq.  2-3. 

3  complex  wave  number  of  the  perturbation  in 

the  z  direction. 

Y  complex  frequency  of  the  perturbation  in 

a  fixed  reference  frame. 

Y'  complex  frequency  of  the  perturbation  in  the 

moving  reference  frame. 

Y*  transformed  complex  frequency  defined  in 

Eq.  2-21. 

9  phase  angle  parameter  defined  in  Eq.  2-12. 

k  amplitude  parameter  defined  in  Eq.  2-19, 

A  angle  of  plane  of  perturbation  with 

respect  to  xy  plane. 

AR  angle  of  resultant  growth  wave  number  vector 


LR 


X_  with  respect  to  x  axis. 


A,  angle  of  oscillation  wave  number  vector  Xy 
with  respect  to  x  axis. 

XR  growth  wave  number  vector. 

Xy  oscillation  wave  number  vector. 

0  wave  number  phase  angle  defined  in  Eq.  2-9. 

0*  wave  number  phase  angle  defined  in  Eq.  2-7. 

ty  wave  number  phase  angle  defined  in  Eq.  2-2. 


I.   THEORETICAL  BACKGROUND  AND  APPROACH 

A.   BACKGROUND 

This  research  deals  with  the  instability  of  plane 
Poiseuille  flow,  that  is,  plane  flow  between  infinite 
parallel  plates.   The  mean  velocity  of  this  flow  is  given 
by  the  expression 


u  =  jd-y2)  .  (l-D 


The  stability  of  such  a  flow  field  is  determined  by 
superimposing  upon  it  an  appropriate  perturbation  and 
determining  whether  this  perturbation  tends  to  grow  or 
decay  over  time.   In  the  present  case  the  perturbations 
are  expressed  by  a  complex  vector  potential  which  is  taken 
to  be  of  the  form 

W  =  [j*G(y)  +  kH(y)]  exp  (ax  +  6z  +  yt)  .       (1-2) 

The  complex  constants  a  and  $  fix  the  spatial  charac- 
teristics of  the  perturbation  and  may  be  arbitrarily  pre- 
scribed whereas  the  complex  constant  y  fixes  the  response  in 
time  and  must  be  found  by  solving  the  vorticity  transport 
equation.  Moreover,  since  a,  6  and  y  are  all  complex,  they 
can  be  resolved  into  real  and  imaginary  components  in  the 
form 
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a  =  aR  +  iaT  (1-3) 


3  =  3R  +  iSj  (1-4) 


Y  =  YR  +  iYj  .  (1-5) 

Thus  the  spatial  characteristics  of  the  perturbation 
are  seen  to  be  completely  defined  by  the  four  constants 
aR,  a,,  3R,  3t  .   In  addition,  the  mean  flow  is  charac- 
terized by  its  Reynolds  number  Re. 

Harrison's  original  analysis  [Ref.  1]  showed  that  the 
perturbation  growth  rate  in  time  as  seen  by  a  fixed  observer, 
and  as  expressed  by  the  parameter  yR,  is  a  definite  function 
of  the  five  parameters  aR,  a,.,  3R,  3T  and  Re,  which  charac- 
terize the  perturbations  and  the  flow.   Thus 


YR  ■  YR[aR,  dj,  3R,  3j,  Re]  .         (1-6) 


B.   PROGRESSIVE  INSTABILITY 

In  a  further  development  of  Harrison's  original  approach, 
Section  II  of  this  thesis  shows  that  three  significant  levels 
of  instability  can  be  defined  which  are  termed  incipient, 
critical  and  fully  developed  instability.   The  definition  of 
these  terms  depends,  in  part,  on  the  algebraic  sign  of  aR. 
However,  Harrison  showed  that  negative  values  of  aR  have  a 
definitely  destabilizing  effect.   Consequently  the  present 
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analysis  is  restricted  to  the  critical  case  of  negative  aR. 
For  this  case  the  three  levels  of  instability  correspond  to 
the  following  levels  of  yR,  namely 

Incipient  Instability  YR)  t  =  0  (1-7) 

Critical  Instability   YR)f  =  -otR  (1-8) 


Fully  Developed  .3  m  en 

Instability  YRJD     I  aR  *      u"yj 


Numerical  solution  of  the  vorticity  transport  equations 
enables  us  to  find  the  three  corresponding  Reynolds  numbers 
at  which  the  above  stability  levels  are  reached.   Thus 

Incipient  Instability  Re)j  =  REj[ctR,  oij,  3R,  Bj]      (1-10) 

Critical  Instability  Re)c  =  Rec[aR,  oij,  3R,  3j]      (1-H) 

Fully  Developed       R  .   =  R   r        R   R  1      ri-121 
Instability  KeJD   ReD^aR'  aI>  BR>  BIJ      U-1ZJ 

C.   TRANSFORMATION  OF  PARAMETERS 

In  Section  II  of  this  thesis,  a  transformation  is 
developed  which  relates  the  original  parameters  aR,  aT,  3R, 


3t,  Re  and  Yd  to  an  alternative  set  of  parameters  which  are 

R' 


*    *        *       * 

symbolized  as  A  ,  0  ,  9 ,  Re   and  y«.   This  transformation 
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can  be  expressed  in  several  alternative  but  equivalent  forms 
In  the  present  context  it  is  convenient  to  write 


aR  =  A*cos(0*+9)  (1-13) 


aT  =  A*sin(0*+9)  (1-14) 


*2 
BR  =  ^j-  {[(l-K2)2+4K2sin26]1/2+cos20*-K2cos2(0*+e)}   (1-15) 

62  =  ^-  {[(l-K2)2  +  4K2sin29]1/2-cos20*+K2cos2(0*  +  6)}   (1-16) 


Re  =  Re*/<  (1-17) 


and 


YR  =  kYr  .  (1-18) 

The  important  fact  about  this  new  set  of  parameters, 
which  are  somewhat  loosely  termed  the  starred  parameters,  is 
that  their  use  permits  the  fundamental  vorticity  transport 
equation  to  be  simplified.   Specifically,  the  relation 
analogous  to  Eq.  1-1  reduces  to 


YR  =  YR[A  ,  0  ,  0,  Re  ]  (1-19) 


The  relations  analogous  to  Eqs .  1-7,  1-8,  and  1-9  reduce 
to 
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Incipient    Instability     Yd)  t    =    °  (1-20) 


Critical    Instability        yR)c   =    "   aR  =    "   f  A  cos C0    +9)         (1-21) 


SJLSrd  VD  ■    -   !  ^  ■    -   I  A*cos(Ae)    (1-22) 


Finally,  the  relations  analogous  to  Eqs .  1-5,  1-6,  and  1-7 
simplify  to 

it  it     A      it 

Incipient  Instability  Re,  =  REy [A  ,0,9]      (1-23) 

it  it     it      it 

Critical  Instability   Rec  =  Rec[A  ,0,6]      (1-24) 


sgbas0^       <  -  **>*>  ?.  ^     ^ 


The  remarkable  feature  of  Eqs.  1-19  through  1-25  is 
that  none  of  these  relations  involve  the  parameter  k. 
Thus  the  number  of  independent  parameters  has  been  reduced 
from  four  in  Eqs.  1-10,  1-11,  and  1-12  to  three  in  Eqs.  1-23, 
1-24,  and  1-25.   This  represents  a  very  significant  simpli- 
fication of  the  problem,  especially  in  view  of  the  tremendous 
computational  burden  which  these  equations  involve. 
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D.   PHYSICAL  SIGNIFICANCE  OF  STABILITY  BOUNDARIES 

The  stability  boundaries  Re  )T,  Re  )f,  and  Re  )D 
symbolized  by  Eqs.  1-23,  1-24,  and  1-25  have  physical 
significance  which  can  be  interpreted  in  a  straightforward 
manner.   Eq.  1-17  shows  that  Re  can  be  regarded  as  the  value 
of  Re  which  corresponds  to  the  reference  case  k  =  1.   Thus 
Re  ) T ,  for  example,  is  the  Reynolds  number  of  incipient 
instability  for  a  perturbation  which  is  characterized  by 

is  ~k 

the  given  values  of  parameters  A  ,  0  ,  and  9  and  by  the 
reference  value  k  =  1.   Since  the  transformed  quantities 
A  ,  0  ,  and  9  may,  at  first,  seem  to  be  somewhat  abstract 
in  character,  it  is  helpful  to  go  back  to  Eqs.  1-13,  1-14, 
1-15,  and  1-16  and  ascertain  the  corresponding  values  of 
the  original  untransformed  parameters  aR,  a-r,  $R,  $t. 

However,  there  is  far  more  to  this  solution  than  just 
the  above  reference  case,  k  =  1.   In  this  connection, 
Eq.  1-17,  when  taken  in  conjunction  with  Eqs.  1-23,  1-24, 
and  1-25,  reveals  a  most  important  result.   It  shows  that 
for  given  values  of  parameters  A  ,  0  ,  and  9,  and  hence 
for  the  corresponding  values  of  Re  )T,  Re  )r,  or  Re  )~, 
the  corresponding  actual  Reynolds  numbers  Re)j,  Re)~,  or 
Re)n  can  be  lowered  indefinitely,  simply  by  increasing 
parameter  k  to  any  desired  extent.   Notice  that  such  a 
shift  of  the  stability  boundaries,  while  it  involves  no 
change  in  parameters  A  ,  0  ,  and  9,  does  involve  changes 
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in  aR,  aT,  3R,  and  3-r.   It  is  therefore  important  to 
summarize  the  nature  of  these  changes  in  the  clearest 
possible  way. 

For  this  purpose,  it  is  useful  to  regroup  the  four 
quantities  aR,  ou,  3R,  and  3y  into  two  vectors  XR  and  Xj 
defined  as  follows: 

XR  =  TaR  +  ¥3R  (1-26) 

I,  =  Tc^  +  Fe.x  (1-2  7) 

(T  and  k  are  unit  vectors  in  the  x  and  z  directions, 
respectively. ) 

Clearly  AR  represents  the  spatial  growth  rate  in  vector 
form,  that  is,  in  terms  of  magnitude  and  direction,  while 
X-r  represents  the  spatial  oscillation  rate  in  like  terms. 
Each  of  these  vectors  is  characterized  by  a  magnitude  and 
a  direction.   In  this  case  the  magnitudes  AR  and  A,,  turn 
out  to  be  governed  by  the  relations 

*? 
7         A  ?  ?    ?    ?   1  /?      7  7    * 

AR  =  ~-  {([(1-k  )z  +  4k  sinze]1/z-(l-K  ))+2cosz0  }       (1-28) 

2A  ????1/?  7  7    * 

Aj   =   *p  {([(1-kz)z  +  4k    sinZ9]i/z-(l-KZ))+2sinZ0    }    .  (1-29) 

Likewise  the  two  corresponding  angles  AR  and  AT  which  the 
above  vectors  make  with  respect  to  the  x  axis  turn  out  to 
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be  governed  by  the  relations 


tan2A  =  [(l-K2)2+4K2sin29]1/2+cos20*-K2cos2(0*+  )      (1.30) 
R  kZ[1  +  cos2(0*+0)] 


*   2.   _  [(l-K2)2  +  4K2sin20]1/2-cos20*  +  K2cos2(0*  +  6)      ri  ,,. 

ta.n     A  j  n 55 \x~ox) 

1  K    [1    -    cos2(0    +0)] 


Thus  the  spatial  form  of  the  perturbations  is  now  fully 
characterized  by  the  four  transformed  parameters  XR,  A,,  AR, 
and  AT  which  are  in  some  respects  more  convenient  than  the 
four  original  parameters  aR,  a,,  8R,  and  3t. 

E.   EFFECT  OF  VARYING  PARAMETERS 

The  connection  between  the  above  perturbation  charac- 
teristics and  the  Reynolds  number  is  still  expressed  by  the 
relation 

Re  =  —  .  (1-32) 

It  is  very  instructive  to  study  the  trends  revealed  by 
Eqs .  1-28  through  1-32  when  parameters  A  ,  0  ,  and  0  are 
held  constant  while  <   is  allowed  to  increase.   Eq.  1-32 
reveals  that  Re)T,  Re)r,  and  Re)n  can  be  decreased  indefi- 
nitely in  this  manner.   On  the  other  hand,  Eq.  1-28  reveals 
that  any  such  decrease  in  Reynolds  number  always  entails  a 
corresponding  increase  in  the  quantity  XR.   Recall  that  AR 
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represents  exponential  growth  rate  in  space.   This  says 
then  that  stability  boundaries  are  not  absolute  in  character 
but  depend  significantly  on  the  "abruptness"  of  the  pertur- 
bation in  space  as  measured  by  parameter  XR.   The  greater 
this  abruptness  parameter,  the  lower  the  Reynolds  number 
at  which  instability  can  occur. 

Conversely,  if  the  permissible  magnitude  of  XR  be 
limited  in  some  definite  manner,  the  reduction  that  can  be 
achieved  in  Re),.,  Re)   and  Re)„  will  be  correspondingly 
limited  as  well.   In  that  case,  a  systematic  exploration 
over  appropriate  ranges  of  the  parameters  A  ,  0  ,  and  9 
should  ultimately  reveal  corresponding  ultimate  stability 
limits  and,  in  particular,  some  critical  Reynolds  number 
below  which  no  instabilities  occur.   Notice,  however,  that 
such  a  critical  Reynolds  number  is  never  absolute,  but  is 
always  contingent  upon  the  restriction  that  has  been  placed 
upon  parameter  XR. 

1.   Restrictions  on  XR 


The  most  obvious  and  direct  restriction  that  can 
be  placed  on  XR  is  simply  to  limit  it  to  some  fixed  value 
or,  for  study  and  comparison  purposes,  to  some  series  of 
successive  fixed  values.   In  general,  the  boundaries  which 
correspond  to  incipient,  critical  and  fully  developed 
instability  will  then  depend  on  the  designated  value  of 
XR.   The  higher  this  value,  the  lower  the  values  of  Re 
at  which  the  above  boundaries  will  occur. 
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Any  such  restriction  of  the  magnitude  of  XR  implies  a 
corresponding  restriction  on  k.   To  show  this,  invert 
Eq.  1-28,  solving  for  k  as  a  function  of  XR.   The  result 
is 


XR\^      2  * 
-  cos  0 


XR\*     .  2.* 

-5r)   +  sin  0 


& 


2    *  2  * 

cos  0   +  sin  0 


(1-33) 


This  relation  may  be  used  in  connection  with  Eq.  1-32 
to  express  the  final  Reynolds  number  at  which  the  designated 
stability  boundary  is  reached.   For  incipient  instability, 
for  example,  this  boundary  may  be  expressed  in  the  form 


Re) 


*   *   * 


_  Re j (A  ,0  ,9)   (1-34) 


4)  +  sin20* 

A 


Analogous  expressions  apply  to  the  boundaries  for  critical 
and  fully  developed  instability. 

Equation  1-34  shows  quite  clearly  that  the  stability 
condition  in  question  depends  on  the  three  characteristic 
parameters  A  ,  0  ,  9  of  the  perturbation  as  well  as  on  the 
limiting  value  assigned  to  parameter  XR.   This  procedure 
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of  calculating  stability  boundaries  for  various  assumed 
combinations  of  A  ,  0  ,  0  and  AR  has  been  carried  out  for 
several  typical  cases  and  the  detailed  results  are  sum- 
marized in  Section  IV  of  this  paper.   Of  course,  these 
examples,  while  representative,  merely  scratch  the  surface 
of  the  stability  problem.   The  complication  remains  that 
the  true  and  ultimate  stability  boundary  represents  the 
lowest  possible  Reynolds  number  at  which  an  instability 
can  just  occur.   This  implies  that  all  possible  combinations 
of  parameters  A  ,  0   and  9  must  be  examined  to  determine 
the  particular  combination  which,  for  a  given  limit  on  XR, 
yields  a  stability  boundary  at  the  lowest  possible  value 
of  Re.   In  other  words,  the  true  stability  boundary  amounts 
to  the  envelope  of  all  the  individual  stability  boundaries. 

Each  individual  boundary  is  characterized  by  some 
specified  combination  of  A  ,  0   and  9  and,  of  course,  also 
of  AR.   Since  there  is  an  unlimited  number  of  such  combi- 
nations, the  amount  of  calculation  involved  in  establishing 
the  desired  stability  envelope  is  prodigious.   Needless  to 
say,  no  such  attempt  was  made  in  the  present  thesis  to 
accomplish  anything  so  ambitious. 

F.   SCOPE  OF  PRESENT  RESEARCH 

The  present  research  was  restricted  to  the  more  modest 

and  realistic  aim  of  calculating  stability  boundaries  for 

ft   ft 
a  few  specific  and  typical  combinations  of  A  ,  0  ,  9  and 

XR.   This  goal  has  been  successfully  attained. 
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G.   PARAMETER  G 

1.   Definition  of  Parameter 

A  detailed  study  of  the  relations  summarized  by 
Eqs.  1-28  through  1-32  reveals  the  possibility  of  express- 
ing a  restriction  on  the  permissible  magnitude  of  AR  in  a 
rather  subtle  and  indirect  way,  through  a  change  of  variable. 
The  particular  algebraic  form  which  the  above  relations 
assume  suggests  the  utility  of  defining  a  new  parameter, 
called  G,  as  follows: 


G2  =  1/2  {[(1-k2)2  +  4K2sin29]1/2  -  (1-k2)}      (1-35) 


This  relation  can  be  readily  inverted  to  give 


K2  -  GW  +  1)  (1.36, 

G   +  sin  9 


2.   Utilization  of  Parameter  G 

Equations  1-35  and  1-36  may  be  used  to  eliminate 
parameter  k  from  Eqs.  1-28  through  1-32,  replacing  it  by 
the  new  parameter  G.   In  this  way  the  following  results  are 
obtained. 

The  vector  amplitudes  XR  and  AT  turn  out  to  be 
related  to  the  new  parameter  G  in  a  fairly  simple  fashion. 
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Thus 

*rfr2       2  *  ,1/2 


XR  =  A  [(G   +  cos"0  )]'  (1-37) 


*    2        7*1/7 
Xj  =  A  [(Gz  +  sin^0  )]1/z  (1-38) 


On  the  other  hand,  the  angles  AR  and  AT  are  not 
simplified  by  the  use  of  parameter  G.   Fortunately  these 
quantities  are  less  significant  than  the  preceding  ones. 
The  governing  equations  become 


tan2A  =  (G2+sin26)(G2-sin20,fe)^G2(G2n)sin2(0%e)      (1_39) 
R  GZ(G2+l)cosZ(0*+9) 


and 


tan2   =  (GZ+sinZ9)(G  -cosZ0  )  +  GZ  (GZ +l)cosZ  (0  +  9)       (1.40) 
1  G  (G  +l)sin  (0  +6) 


The  important  Reynolds  number  relation  below  is  once 
again  simple.   For  definiteness  it  is  written  specifically 
for  the  case  of  incipient  instability,  by  analogy  with 
Eq.  1-34.   Similar  expressions  apply  also  to  the  boundaries 
of  critical  and  fully  developed  instability.   Thus 
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Re)   =  / %  V   9  ReTCA  ,0  ,9)  (1-41) 

JGZ(GZ  +  1)    X 


Equation  1-41  shows  that  the  stability  depends  on 
the  particular  parameters  A  ,  0  ,  9  and  on  the  limiting 
value  assigned  to  G.   Hence  G  plays  a  similar  role  in 
relation  to  Eq.  1-41  that  XR  plays  in  relation  to  Eq.  1-34. 

The  results  summarized  elsewhere  in  this  thesis 
are  presented  primarily  from  the  perspective  expressed  by 
Eq.  1-34.   The  alternative  version  shown  by  Eq.  1-41  is 
included  in  this  discussion  because  of  its  theoretical 
interest,  but  this  version  is  not  used  in  the  presentation 
of  calculated  results. 

Notice  that  in  either  version,  assuming  some  assigned 
limit  for  XR  or  G,  an  exploration  is  still  required  over  the 
domain  of  parameters  A  ,  0   and  9  to  find  the  particular 
combination  that  yields  the  minimum  Reynolds  number.   Of 
course,  such  extensive  exploration  could  not  be  undertaken 
in  the  present  thesis  owing  to  time  limitations. 

H.   REDUCTION  TO  CLASSICAL  THEORY 

It  is  pertinent  to  note  that  the  classical  theory  of 
the  stability  of  plane  Poiseuille  flow  amounts  to  a  special 
case  of  the  more  general  theory  discussed  above.   It 
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amounts,  in  fact,  to  the  special  case  for  which 

AR  =  0.  (1-42) 

Study  of  Eq.  1-28  reveals  that  Eq.  1-42  can  be  satisfied 
if  and  only  if  we  set 

9=0.  (1-43) 

and 

0*  =  \   .  (1-44) 

From  Eqs .  1-28,  1-32,  and  1-42  we  may  infer  also  that 

k  =  1  .  (1-45) 

It  then  follows  from  Eq.  1-35  that 

G  =  0. 

Moreover,  we  also  find  under  these  conditions  that 

aR  =  0  ,  3R  =  0  ,  and  6X  =  0.  (1-46) 

It  is  evident  that  the  general  theory  discussed  in  this 
thesis  is  immensely  more  comprehensive  than  the  classical 
theory  as  limited  by  Eqs.  1-42  through  1-46. 
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II.   SIMILARITY  TRANSFORMATION 

The  perturbation  characteristics  are  fully  defined,  as 
in  Ref.  1,  by  the  four  real  parameters  aR,  a,,  3R,  and  3j. 
ctR  and  cu  are  the  components  of  the  complex  wave  number  of 
the  perturbation  in  the  x  direction  and  3R  and  3,  are  the 
components  of  the  complex  wave  number  of  the  perturbation 
in  the  z  direction. 

The  above  parameters  satisfy  the  following  relations: 


*  10 
a  =  kA  e    =  aR  +  iaT  (2-1) 


3  =  oA*ellp   =  3R  +  iSj  (2-2) 


A  very  useful  alternative  set  of  parameters  is  ctR,  a,,  9 
and  k. 

*2 
a   is  defined  by 


*2     7 
a   »  oT  +  3   .  (2-3) 


*    *    *      * 

A  ,  0  ,  aR  and  aT  are  defined  by 


* 

a   =  A  e  v      =  aR  +  ia-j-  .  (2-4) 
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aR  and  aT  are  the  components  of  a  ,  the  transformed  complex 
wave  number  parameter  and  k  is  the  transformed  perturbation 
amplitude  parameter. 

k  and  9  are  defined  by 


*  ii 

a  =  a  Ke 


(2-5) 


A  and  <   are  positive  by  definition. 

Given  the  original  parameters  aR,  a,,  $R  and  3T,  the 
transformed  parameters  aR,  a-.,  0  and  k  may  be  deduced  from 
equations  2-1  through  2-5  as  follows: 


A*  =  [(a2  -  a2  +  32  -  62)   +  4(aRaI  +  3R3j)2]        (2-6) 


0   =  1/2  arctan 


2^aRaI  +  W 


aR  *  aI  +  3R  "  3I 


(2-7) 


kA*  =  [a2  +  a2]1/2 


(2-8) 


/a 
0  =  arctan/  — 
\aR 


(2-9) 


*     *      * 

aR  =  A  cos  0 


(2-10) 
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X        A  A 

a,  -  A  sin  0  (2-11) 


=  (0  -  0*)  (2-12) 


The  governing  equations  of  Ref.  1,  using  the  five 
independent  parameters,  Re,  aR,  ou ,  3R,  3-r  are  as  follows: 
Using  the  auxiliary  expressions  below, 


2     2 
T(7)  =  a  R*  B   -  a|  (1-y2)  (2-13) 


2  + 

"Ee~ 


Uy)  =  a  n!  3   +  T(y)  (2-14) 


J(y)  =  (a2  +  $2)T(y)  -  3a  (2-15) 


the  fundamental  vorticity  equation  becomes 


[Jg.  Hlv  +  I(y)H"  ♦  J(y)H]  -  Y[H"  +  (a2  +  32)H]  =  0.    (2-16) 


The  associated  vorticity  equation  is 


[Re  G"  +  (T(y)-Y)G  =  -^r-T[LHM'  +(T(y)-Y)H'-3ayH]  .  (2-17) 

a  +3 


The  number  of  independent  parameters  in  equations  2-13 
through  2-17  can  be  reduced  to  four  by  utilizing  the 
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following  functional  transformations 


Re  =  K-1Re*  (2-18) 


2     2*^ 
a      +   B   =  a  (2-19) 


*  i8 

a  =  a  Ke  (2-20) 


Y  =  <Y*  (2-21) 


H(y)  =  K~Vl9H*(y)  (2-22) 


G(y)  =  ic"1e"ie0G*Cy]  (2-23) 


Substitution  of  Eqs.  2-18  through  2-23  into  the  general 
solution,  Eqs.  2-13  through  2-17  yields  the  three  auxiliary 
functions 

*2 
*       a,       *  i8  ^     2 
T  (y)  =  %-  -  a  e1"  i(l-y^)  (2-24) 

Re  L 


*2 

I"(y)  =  ^-r  +  T"(y)  (2-25) 

Re 


*      a      * 


*2 


J  (y)  =  a   T  (y)  -  3a  e   .  (2-26) 
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The  principal  vorticity  transport  equation  becomes 

[i_^H    +1  (y)H   +J  Cy)H  ]  -  Y  [H   +a   H  ]  =  0  (2-27) 

Re 

The  associated  vorticity  transport  equation  becomes 


1       *M      *  *        1       1     *!tt       *  *     *t 

[i-TfC   +(T  (y)-Y)G   =  -^t-Ai«    +(T  (y)-Y  )H 

Re  *L   Re 

a 

*  ie   * 
-3a  e  yH  ]  .      (2-28 


Equations  2-24  through  2-28  now  involve  only  four 

•k  *     * 

parameters  Re  ,  aR,  a-r  and  0.   k,  the  fifth  parameter,  has 
cancelled  out.   k  becomes  part  of  the  solution  again  during 
the  reverse  transformation  of  results  from  starred  parameters 
to  the  original  parameters. 
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III.   STABILITY  CRITERION 

A.  BACKGROUND 

Studies  of  the  stability  of  Poiseuille  flow  have  used 
various  criteria  for  determining  the  stability  of  the  flow 
from  the  solutions  obtained.   The  growth  rate  in  time,  yR» 
is  usually  used  when  there  is  no  real-exponential  spatial 
variation  [Salven  and  Grosch,  1972].   When  exponential 
growth  in  space  has  been  included  [Garg  and  Rouleau,  1971] 
the  real  part  of  the  spatial  wave  number  has  been  used  to 
give  the  instability  but  this  procedure,  while  seemingly 
plausible  at  first  inspection,  cannot  be  really  justified 
with  any  rigor.   The  Lagrangian  approach  described  below 
is  believed  to  be  a  superior  method  for  dealing  with  this 
case.   In  other  cases,  stability  has  arbitrarily  been 
evaluated  with  respect  to  a  frame  of  reference  moving  down- 
stream at  the  phase  velocity  of  the  perturbation  but  again, 
this  procedure  has  no  strict  rational  justification,  and 
especially  so  in  connection  with  the  fully  three-dimensional 
perturbations  considered  in  this  thesis. 

B.  LAGRANGIAN  REFERENCE 

For  perturbations  that  are  both  oscillatory  and  have 
exponential  rates  of  growth  or  decay  in  the  streamwise  and 
transverse  directions  a  Lagrangian  frame  of  reference 
proves  useful.   The  fluid  particles  have  velocities  varying 
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from  0  to  1.5  depending  on  their  distance,  y,  from  the 
walls.   The  velocity  distribution  is  given  by 


U  =  f(l-y2)  .  (3-1) 


Consider  a  coordinate  system  moving  in  the  x  direction 
with  the  mean  velocity  of  a  given  fluid  particle.   Let  y 
be  the  mean  vertical  coordinate  of  the  moving  particle. 
Then  the  velocity  of  the  moving  reference  frame  is  the  same 
as  the  velocity  of  the  streamline  along  which  the  above 
particle  moves  and  is  given  by  Eq.  3-1.   Let  x',  y,  z,  t 
be  the  coordinates  and  a,  3,  y'  the  complex  wave  numbers  with 
respect  to  the  moxing  axes.   The  form  of  the  perturbation 
vector  potential  for  a  given  eigenvalue  obtained  as  a 
solution  is 

W  =  TG(y)  +  k~H(y)  exp  (ax  +  $z  +  yt)  .        (3-2) 

The  complex  frequency  y'  seen  from  this  moving  reference 
is  different  than  from  a  fixed  frame.   To  relate  y '  to  y  the 
perturbation  vector  potential  is  written  in  the  moving 
frame  and  transformed  into  the  form  for  the  fixed  frame. 
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W  =  (jG  +  kH)exp(ax'  +  Sz  +  y't) 

=  (jG  +  kH)exp(a(x-Ut)  +  gz  +  y't) 

=  (jG  +  kH)exp(ax  +  6z  +  (y '  -  aU)t) 

=  (j~G  +  kH)exp(ax  +  6z  +  yt)  (3-3) 

Therefore 

yf  -  aU  =  y  (3-4) 

Solving  for  y'  and  splitting  into  real  and  imaginary 
parts  yields 

YR  =  YR  +  aRU  (3"5) 


YJ-  =  Yj  +  ctjU  (3-6) 


If  yX   is  positive,  zero,  or  negative,  the  perturbation  is 
said  to  be  unstable,  neutral,  or  stable,  respectively,  with 
respect  to  the  moving  reference  frame.   Thus,  the  value  of 
y'  is  taken  to  be  a  measure  of  the  stability  of  each  eigen- 
value obtained. 
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C.   TRANSFORMED  STABILITY  CRITERION 

Consider  now  the  transformation  to  the  starred  parameters 
The  condition  of  stability  is  determined  by  the  value  of  y' 
which  Eq.  3-5  gives  as 

YR  =  YR  +  aRU  '  (3'5) 


Now 


t 

YR  =  KYR 


* 

YR  =  <YR  (3-7) 


aR  =  <aR 


Substitution  of  Eqs.  3-7  into  Eq.  3-5  yields 


YR   =  YR  +  aRU  (3-8) 


aR  is  defined  by  Eq.  2-10  as 


•*•     *      a 
aR  =  A  cos(0  +9)  (2-10) 


Therefore 


*     *    *     * 


YR  =  YR  +  A  cos(0  +9)  (3-9) 
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The  three  stability  boundaries,  incipient,  critical, 
and  fully  developed  are  determined  by  the  condition  that 

exists  when  yR  =  0.   Setting  Eq.  3-9  equal  to  zero  and 

* 

solving  for  yD  yields 


R 


YR  =  UA*cos(0*+9)  .  (3-10) 


Now  incipient  instability  corresponds  to  zero  growth  rate 
with  respect  to  a  coordinate  system  which  moves  with  the 
flow  velocity  at  the  wall,  which  is  zero. 


y*R)1    =  0  (3-11) 


Critical  instability  corresponds  to  zero  growth  rate  with 
respect  to  a  coordinate  system  which  moves  with  the  mean 
velocity  of  the  flow,  which  is  unity. 


YR)C  =  -A  cos(0  +9)  (3-12) 


Fully  developed  instability  corresponds  to  zero  growth 
rate  with  respect  to  a  coordinate  system  which  moves  with 
the  velocity  of  the  flow  on  the  centerline. 


Yr)d  =  -  jA*cos(0*+0)  (3-13) 
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IV.   RESULTS 

A.   TRANSFORMED  PARAMETERS 

Equations  2-24  through  2-28  were  solved  on  the  IBM-360 

R 


* 

computer  and  the  most  unstable  growth  rate,  yR,  for  a  given 


0  ,  9 ,  Re  ,  A  was  obtained.   The  boundaries  for  incipient 
and  critical  instability  were  determined  by  the  criteria 
explained  in  Section  II.   Fully  developed  instability  did 
not  occur  for  the  case  studies.   Two  values  of  0  were 
studied  at  various  values  of  8,  Re  ,  and  A  . 
1.   0*  =  90° 


Values  of  9  explored  were  0,  1,  2,  3,  4,  5,  and  6°. 
Because  one-degree  increments  of  9  yield  graphical  results 
that  are  extremely  cluttered,  this  study  will  present  the 
results  only  for  9  =  0,  3,  and  6°.   Figure  4-1  shows  the 
stability  boundaries  for  0   =  90°.   The  effect  of  changing 
9,  while  holding  0   constant,  can  be  seen.   An  increase  in 
9  causes  a  degrease  in  Re   for  both  incipient  and  critical 
instability.   Also,  for  a  given  9,  the  boundary  for  critical 
instability  occurs  at  a  higher  Re   than  does  the  boundary 
for  incipient  instability. 

There  is  only  one  curve  for  9=0°.  In  this  case 
the  criteria  for  incipient,  critical,  and  fully  developed 
instability  turn  out  to  be  identical.   Equation  3-10  shows 


Yp  =  -UA*cos(0*  +  9)  (3-10) 
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and  cos  (90°  +  0°)  =  0;  thus  the  criterion  for  stability  on 
the  three  boundaries  is  yR  =  0. 

Note  that  the  boundary  of  incipient  instability  for 
0=3°  shows  an  abrupt  distoncinuity  represented  by  segment 
ab.   This  discontinuity  is  similar  to  that  obtained  by 
Harrison.   It  is  expected  that  extension  of  the  incipient 
boundary  for  other  values  of  9  would  reveal  the  same 
characteristic. 

Figures  4-2  and  4-3  are  plots  of  the  growth  rate, 
YR>  versus  Re   for  9=3°  and  6  ,  respectively.   Both  plots 
show  the  locus  of  points  that  represent  the  boundaries  of 
incipient  and  critical  instability.   It  can  be  seen  that 
fully  developed  instability  is  not  reached  even  at 
Re*  =  100,000. 


2.   0   =  95 


o 


Due  to  a  lack  of  time  only  two  values  of  9  were 
explored,  9=0  and  3°.   A  comparison  of  Fig.  4-4  with  4-1 
shows  that  the  stability  contours  follow  much  the  same 
pattern  for  both  cases.   However,  increasing  0   to  95° 
causes  a  corresponding  decrease  in  Re  . 

Figures  4-5  and  4-6  show  the  growth  rate  as  a 
function  of  Re   for  0=0  and  3°,  respectively.   The  locus 
of  points  that  represent  the  boundaries  of  incipient  and 
critical  instability  can  again  be  seen.   Fully  developed 
instability  is  not  reached  at  Re   =  100,000. 
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B.   RESULTS  TRANSFORMED  TO  UNSTARRED  PARAMETERS 

In  order  to  present  the  results  in  a  manner  consistent 
with  Ref.  1,  it  is  necessary  to  transform  the  results  from 
starred  to  unstarred  parameters.  This  transformation  puts 
the  results  in  a  more  easily  understandable  form. 

From  Fig.  4-7  the  following  relations  can  be  deduced. 


-jr]  =  1/2(k2  +  [(1-k2)2  ♦  4K2sin29]1/2+  cos20*)     (4-1) 


>x  v2 

i)    =    1/2(k2    +    [(1-k2)2    +    4<2sin26]1/2-    cos20*)  (4-2) 


tan2A      =    [U-<V  +  4K2sin29]1/2  +  cos20*-K2cos2(0*  +  9)  (4_3) 

R  k2[1    +   cos20X+0) ] 


72  2  21/7  *?  * 

tan2A      =    [(!-<)    +4<    sin   9]1/^-cos20    +<   cos2(0    +9) 

1  <Z[1    -    cos2(0X+9)] 


(4-4) 


Aj.  and  AR  are  the  magnitude  and  direction,  respectively, 
of  the  perturbation  and  growth  vector.  X-r   and  AT  are  the 
magnitude  and  direction,  respectively,  of  the  perturbation 
oscillation  vector. 

1.   Perturbation  Rate  Vectors ,  Magnitude 

Figure  4-8  shows  G  as  a  function  of  Re/Re  with  9 
as  an  independent  parameter;  the  curves  are  valid  for  all 
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Figure  4.7.   Vector  Diagram  for  Parameter  Trans formati 
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values  of  0  .   For  any  fixed  value  of  0,  G  decreases  as 
the  Reynolds  number  ratio  increases.   In  the  case  of  9  =  0° 
G  decreases  with  increasing  Reynolds  number  ratio  until 
Re/Re   =1.   G  then  maintains  a  constant  value  of  zero  for 
further  increases  in  the  Reynolds  number  ratio.   Although 
values  of  9  above  6°  were  not  explored,  values  of  9  up  to 
60°  are  shown  here  to  demonstrate  the  trend  as  9  increases. 

2.   Perturbation  Growth  Rate  Vectors,  Direction 

2  2 

The  quantities  of  tan  AR  and  tan  AT  are  fixed  by 

2 
Eqs .  4-2  and  4-3,  respectively.   However,  if  tan  AR  be 

specified,  this  does  not  fix  AR  uniquely  as  there  are 

four  angles,  one  in  each  quadrant,  which  have  the  specified 

value  of  tan  AR.   Similar  considerations  apply  also  to  the 

other  angle  A,.   Hence,  to  determine  AR  and  A.  uniquely 

it  is  necessary  to  consult  auxiliary  relations  which  fix 

the  quadrant  in  which  these  angles  really  fall. 

The  components  of  AR  which  fix  angle  AR  are 


aR  =  XRCOsAR  =  kA  cos(0  +e)  (4-5) 


and 


The 


3R  =    XRsinAR   =   aA  cosijj    .  (4-6) 


components  of  A,  which  fix  angle  A,  are 


a j  =  AjCOsA,  =  kA  sin(0  +9)  (4-7) 
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and 


* 


3j  =  AjSinAj  =  aA  sinij;  .  (4-8) 


Moreover,  in  this  study,  the  angle  (0  +9)  has  been 
restricted  to  lie  in  the  second  quadrant  so  that 


aR  <  0  (4-9) 


a j  >  0  (4-10) 

Recall  that  negative  values  of  aR  were  shown  by 
Harrison  to  be  destabilizing.   That  is  why  the  present 
study  is  restricted  to  negative  values  of  aR. 

In  order  to  determine  the  algebraic  signs  of 
components  8R  and  Br,  it  is  necessary  to  bracket  the  range 
of  the  angle  ty .      A  study  of  Fig.  4-7  reveals  that  for 
positive  values  of  9,  the  following  limits  apply. 


lim  =  0*  (4-11) 

-►  0 


* 


IT 


lim  =  (0  +9)  -  y  (4-12) 

k  -*  0  z 


It  is  also  evident  that  the  x-y  plane  is  a  plane  of 
symmetry  and  that  therefore  a  reversal  of  the  perturbation 
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characteristics  with  respect  to  the  z  axis  is  permissible 
and  leaves  the  essential  features  of  the  solution  other- 
wise unchanged.   This  amounts  to  saying  that  the  angle  ty 
can  always  be  changed  by  ±180°,  with  no  significant  effect 
upon  the  solution  except  for  a  reversal  of  the  perturbations 
with  respect  to  the  plane  of  symmetry.   For  definiteness 
in  this  discussion,  however,  we  limit  the  angle  ty    as  inci- 
cated  by  Eqs .  4-9  and  4-10.   It  is  then  evident  that  the 
above  auxiliary  relations,  along  with  the  basic  relations 

of  Eqs.  4-3  and  4-4  ,  now  suffice  to  fix  AR  and  AT  uniquely 

* 
for  any  assigned  values  of  the  parameters  0  ,  8  and  k. 

Figure  4-9  shows  AR  as  a  function  of  Re/Re   for 
0   =  90°  with  9  as  an  independent  parameter.   This  plot 
shows  a  constant  value  of  AR  ■  90°,  for  6=0°  and  Re  <  Re  . 
When  Re  >  Re  ,  the  vector  magnitude,  AR  is  zero.   For  all 
other  values  of  9  the  perturbation  growth  rate  vector,  AR, 
rotates  from  near  the  transverse  to  near  the  upstream  direc- 
tion  as  Re/Re   increases. 

Figure  4-10  shows  the  rotation  of  the  perturbation 
growth  rate  vector  as  Re/Re   increases,  for  0   =  95°.   Note 
that  when  0   =  90°,  AR  varies  between  90°  and  180°  whereas 
when  0   =  90°,  AR  varies  between  90°  and  270°. 

3.   Perturbation  Oscillation  Rate  Vectors,  Direction 

Figures  4-11  and  4-12  are  similar  plots  showing 
the  rotation  of  the  oscillation  rate  vector  with  changing 
Re/Re   for  0   =  90°  and  95°,  respectively.   Comparing 
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Fig.  4-11  with  4-9,  for  0   =  90°  and  0=0°  both  AR  and  Aj 
have  constant  values  when  Re/Re   <  unity. 

4.   Stability  Boundaries  in  Unstarred  Parameters 

Figures  4-13  and  4-14  show  the  stability  boundaries 
for  unstarred  parameters.   A   is  0.05  for  both  plots.   A 
comparison  with  Figs.  4-1  and  4-4  shows  that  the  character 
of  the  boundaries  remains  unchanged.   However,  Reynolds 
number  is  greater  than  starred  Reynolds  number. 

Although  it  is  not  shown  here,  it  was  found  that 
increasing  AR  causes  the  boundaries  to  move  to  the  left 
and  up.   In  other  words,  an  increase  in  AR  causes  an 
increase  in  A,  and  a  decrease  in  Reynolds  number. 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  research  reported  here  is  based  largely  on  the 
theory  developed  by  Harrison  in  Ref.  1.   However,  this 
theory  is  further  extended  in  the  present  work  by  the 
introduction  of  a  useful  similarity  transformation.   The 
transformation  reduces  the  number  of  independent  parameters 
required  in  the  computer  solution  from  five  to  four  and 
thereby  significantly  reduces  computation  time. 

In  Ref.  1,  Harrison  showed  that  negative  values  of  aR 
are  destabilizing.   The  results  presented  here  show  that, 
for  negative  values  of  aR,  the  critical  Reynolds  number  for 
plane  Poiseuille  flow  can  be  lowered  indefinitely  by 
proper  selection  of  perturbation  characteristics,  provided 
that  the  corresponding  increase  in  AR  is  acceptable.   Thus 
the  stability  boundaries  are  not  absolute  in  character  but 
depend  significantly  on  the  "abruptness"  of  the  perturbation 
in  space  as  measured  by  parameter  AR.   The  lowering  of  the 
critical  Reynolds  number  in  this  way  provides  one  possible 
explanation  for  the  disagreement  between  earlier  theory 
and  experiment. 

The  stability  of  the  flow  along  a  particular  streamline 
has  been  shown  to  depend  on  its  velocity.   Negative  values 
of  aR  yield  the  greater  instabilities  and  streamlines  with 
the  lowest  velocities,  those  nearest  the  walls,  will  be 
the  most  unstable.   This  seems  to  agree  with  experiment. 
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According  to  Schlichting,  transition  from  laminar  to 
turbulent  flow  is  "characterized  by  an  amplification  of 
the  initial  disturbances  and  by  the  appearance  of  self- 
sustaining  flashes  which  emanate  from  fluid  layers  near 
the  wall  along  the  tube." 

Instability  was  found  to  be  progressive  in  nature  and 
two  of  the  three  defined  stability  boundaries  were  located, 
incipient  and  critical.   Further  research  needs  to  be  done 
for  a  wider  range  of  parameters  A  ,  0  ,  and  9  to  find  the 
combinations  that  correspond  to  minimum  Reynolds  numbers  at 
the  above  stability  levels. 
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APPENDIX  A 
USE  OF  THE  COMPUTER  PROGRAM 

It  was  found  extremely  useful  to  precompile  the  program 
on  a  disk  thereby  avoiding  the  inconveniences  of  reading  in 
the  complete  card  deck  for  each  run.   An  additional  advantage 
is  a  reduction  in  turn-around  time  of  considerable  magnitude. 
The  following  will  give  procedures  and  hints  that  can  be 
found  in  the  W.  R.  Church  Computer  Center  but  require  time- 
consuming  search. 

1.  Pre-compiling  Program 

To  compile  the  program,  the  following  was  read 
into  the  system 
//  Green  Job  Card, Time  =(0,59) 
//  EXEC  FORTCL 
//  FORT.SYSIN  DD  * 

/* 

Program  Card  Deck  goes  here.   (no  data) 
//LINK.SYSLMOD  DD  DSNAME-S2593 . LIB (POIS) ,DISP= (NEW, KEEP) 
//  UNIT=2  321,VOLUME=SER=CEL006,LABEL=RETPD=2  20, 
//  SPACE=(CYL, (6,1,1) ,RLSE) 

/* 

2.  Program  Execution 

Once  the  program  is  compiled,  running  the  program 
consists  of  punching  data  cards  in  the  namelist  format  and 
reading  them  in  with  the  following  deck  of  cards. 
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//  Green  Job  Card,Time= (0 , 59) 

//  GO  EXEC  PGM=POIS,REGION=178K 

//STEPLIB  DD  DSNAME=S2593.LIB,DISP=SHR, 

//  VOLUME=SER=CEL006,UNIT=2321 

//FT06F001  DD  sysout=A,DCB= (RECFM=FBA,LRECL=133,blksize=3325) 

//FTOSF001  DD  * 

SLIST  N=30,REY=3000,TH=.0  5  2360,ASTAR=1.8,PHIS=95,$END 

/* 

Note:   Column  1  is  blank  on  the  list  card. 

Three  decks  of  these  cards  were  used  with  each  having  a 
different  job  name,  i.e.,  NEWBY  64A,  NEWBY  64B,  and  NEWBY  64C 
This  proved  useful  as  three  jobs  could  be  loaded  at  one  time 
and  one  could  keep  track  of  what  was  already  printed  and 
what  remained  to  be  processed.   It  was  also  found  to  be 
useful  to  have  three  sets  of  job  cards  with  each  set  having 
a  different  time.   For  Time=(0,59),  59  sec,  one  list  card 
(data)  was  inserted.   This  was  used  for  quick  turn-around 
time  and  only  a  few  points  were  being  explored.   For 
Time=(2,00),  2  minutes,  three  list  cards  could  be  read  in 
and  for  Time=(4,00)  six  list  cards  could  be  used.   Occasion- 
ally the  four-minute  time  parameter  would  terminate  execu- 
tion after  five  list  cards  had  been  processed. 

3.   Program  Alteration  after  Compilation 

If  changes  were  to  be  made  in  the  program  the  file 
was  scratched  and  a  new  file  established  with  the  changes 
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incorporated.   To  scratch  the  program  on  file  the  following 
deck  was  used. 

//  Green  Job  Card 

//  EXEC  PGM=IEHPROGM 

//SYSPRINT  DD  SYSOUT=A 

//DD1  DD  UNIT=2321,VOL=SER=CEL006,DISP=SHR 

//SYSIN  DD  * 

SCRATCH  VOL  =  2  321  =  CEL006,DSNAiME  =  S2  59  3.LIB,PURGE 
/* 

Note:   The  scratch  card  begins  in  column  5. 

In  all  three  decks  the  name  of  the  program  (POIS) 
appears.   The  choice  of  a  program  name  is  an  individual 
choice  but  once  chosen  it  must  appear  the  same  in  all 
card  decks.   The  only  other  item  that  appears  with  unique- 
ness is  the  individual  user  number.   In  the  context  of  this 
paper  that  number  was  2593  and  must  agree  with  the  user 
number  on  the  job  card. 

There  are  two  possible  selections  on  input  parameters 
for  obtaining  data.   Both  are  used  in  this  study.   It  is 
possible  to  select  values  of  0  ,  0,  Re  ,  and  vary  A   to 
find  a  point.   This  method  was  used  first  and  worked  well 
when  obtaining  a  solution  from  the  computer.   The  problem 
arises  when  interpreting  and  transforming  the  results.   It 
proves  useful  to  have  values  for  fixed  A  and  this  method 
does  not  provide  this  easily. 
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An  alternative  technique  is  to  select  0  ,  6,  A  ,  and 
then  vary  Re  .  To  construct  Figures  3-1  through  3-6  the 
fixed  A   technique  provides  data  in  an  easily  usable  form 
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APPENDIX  B 
CHANGES  TO  COMPUTER  PROGRAM  IN  REFERENCE  1 

The  following  changes  were  made  to  the  computer  program 
in  Ref.  1  to  convert  to  starred  parameters. 

1.  Program  #1 

Statement  number  0002  (COMPLEX  *16  A,B)  was  deleted 
and  two  type  declaration  statements  (REAL*8  TH)  and 
(C0MPLEX*16  A)  were  inserted.   The  namelist  statement, 
number  0008,  was  revised  to  read:   NAMELIST  /  LIST  /N,REY, 
TH,ASTAR,PHIS,VEL.   Statement  number  0018,  B  =  DCMPLX(BR,BI) 
was  deleted  and  the  following  statements  added:   PHI  = 
PHIS/57.2958,  AR  =  ASTAR  *  (DCOS(PHI)),  AI  =  ASTAR  *(DSIN 
(PHI)),  THD  =  (TH*180.0)/3. 141592654.    Other  changes  to 
program  #1  were  those  required  to  write  out  the  revised 
inputs . 

2.  Subroutine  DEIGEO 

One  small  change  was  made  to  this  subroutine  due  to 
the  fact  that  values  required  by  external  functions  CHM1E1 
and  CHM2E1  were  passed  by  DEIGEO.   Statement  0010  of  DEIGEO, 
B  =  BETA,  was  deleted  and  TH  =  THETA  was  inserted.   The 
common  statement  and  type  declaration  statements  were 
revised  to  incorporate  the  change  to  starred  parameters. 

3.  Functions  CHM1E1  and  CHM2E1 

External  functions  CHM1E1  and  CHM2E1  required 
extensive  modification  as  follows: 
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The  type  declaration  statement  (REAL*8  TH,DUR)  was  added, 
CH4M1(Y)  =  A/REY  was  changed  to  CH4M1 (Y)  =  A*EI/REY. 
CH2Ml(y)  was  changed  to  equal 

-1.5DO*A**2*EI2*(lDO-y**2)+2DO*AEI*(A**2)/REY  . 
CHOMl(Y)  was  changed  to  equal 

-AEI*((A**2)*(1.5D0*AEI*(1D0-Y**2)-(A**2)/REY) 
+3D0*AEI) 
CH2M2(Y)  =  A  changed  to  CH2M2(Y)  =  AEI. 
The  following  statements  were  added  after  CH2M2(Y)  and 
ENTRY  CHM2El(k,Y) : 

DUR  =0.0 

DU  =  DCMPLX(DUR,TH) 

EI  =  CDEXP(DU) 

AEI  =  A*EI 

E12  =  CDEXP(2*DU) . 
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PROGRAM  ¥1 


PROGRAM  TO  PRINT  EIGENVALUES 
FOR  THE  3-0  POISEULLE  FLCW  PROBLEM 


THIS  PROGRAM  SOLVES  THE  LINEARIZEC  NAVIER  STOKES 
ECUATIQN  FOR  POISEULLE  FLOW.   THE  EIGENVALUES 
RESULTING  FROM  THE  FINITE  DIFFERENCE  APPROXIMATION 
AFE  PRINTED. 


INPUT 


THE  FOLLOWING  MAY  BE  INPUT  TO  THE  PROGRAM  *S  DATA 
USING  NAMELIST,  'LIST'.   NOTE,  THE  DEFAULT  VALUES 
ARE  ONLY  SET  INITIALLY  AND  VALUES  SET  BY  THE 
USER  WILL  NOT  BE  CHANGED  BETWEEN  RUNS 

N  -  HALF  OF  THE  NUMBER  OF  FINITE  DIFFERENCE  GRID 
POINTS  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  END 
POINTS.   N  MUST  BE  .LE.  MDIMt  WHICH  IS  THE 
DIMENSION  OF  THE  MATRICES  IN  THIS  PROGRAM. 
DEFAULTED  TO  THE  VALUE  OF  NDIM,  THAT  IS  THE 
DIMENSION  OF  THE  MATRICES.   SEE  PROGRAM  BELOW 
THE  DEFAULT  VALUE. 


FOR 


REY  -  THE  *  REYNOLDS  NUMBER 
VALUE  =  60C0.0 


(<?EAl*8)  DEFAULT 


AR, AI  -  THE  REAL 
THE  STARRED  WAVE 
DEFAULTED  TO  0.0 


AND  IMAGINARY  PARTS 
NUMBERS   <3EAL*8) 
AND  1.0  RESPECTIVELY 


OF 


VEL  -  THE  VELOCITY  OF  THE  MOVING  COORDINATE 
REFERENCE  SYSTEM  FOR  WHICH  THE  STABILITY  IS 
DETERMINED.   (REAL*8)   DEFAULTED  TO  3.0 


OUTPUT 


THE  OUTPUT  OF 
EIGENVALUES. 


THIS  PROGRAM  IS  A  TABULATION  OF  THE 
TWO  LISTS  ARE  PRINTED,  ONE  FOR  THE 
EIGENVALUES  CORRESPONDING  TO  EVEN  E I  GEN FUNCT I GNS 
AND  ONE  FOR  THOSE  CORRESPONDING  TO  DDD  EIGEN- 
FUNCTIONS.   THE  STABILITY  OF  EACH  EIGENVALUE  IS 
PRINTED  AND  THE  LEAST  STABLE  EIGENVALUE  IS  MARKED 
WITH  ASTERISKS.   A  PLOT  OF  THE  EIGENVALUES  IS  ALS^ 
PRINTED. 

SUEROUTINES 

THIS  PROGRAM  CALLS  THE  SUBROUTINE  •DEIGEO'  TO 
SOLVE  FOR  THE  EIGENVALUES.   SUBROUTINE  'PLCTP' 
IS  USED  TO  PLOT  THE  EIGENVALUES  <2H    THE  PRINTER. 


IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  TH 

CCMPLEX*16  A 

REAL*4  GR4(60) ,GI4(60I 

REAL*3  GRE<60) ,GIE(60) ,GR0(60) , GIC(oO) 

C0MPLEX*16  XMAT<60,30,3) 

CO  MP  LEX* 16  YM AT (60, 30) ,WVEC(60) ,BMAT( 5,60) 
EQUIVALENCE ( Y MAT (1, 1),XMAT(1,1,3)), 

*  ( SMAT(l,l),XMAT(i, 1,3) ) , 

*  (WVEC(l)  ,XMAT( 1,6,3) ) 
NAMELIST  /  LIST  /  ,N  ,  REY  ,TH,  AST  A  R  ,  PHI  S  ,  VE  L 
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C  INITIALIZE  VARIABLES   (SET  DEFAULT  VALUES) 
C 

MOIM  =  60 

N  =  60 

REY  =  6000D0 

Th  =  000 

VEL  =  0C0 
C 

C  READ  NAMELIST  AND  SET  ALPHA  AND  ASTAR 
C 

1  REA0(5,LIST,ENO=100) 

PHI  =  PHIS/57.2958 

AR  =  ASTAR  *(DCOS(PHI )) 

AI  =  ASTAR  *  CDSINCPH1 ) ) 

A  =  DCMPLX(AR,AI ) 

TFD  =  (TH*180.0)/3. 141592654 
C 

C  PRINT  INPUT  VALUES  AS  PAGE  HEADI.MG  FOR  EIGENVALUE  LIST 
C 

WRITE<6,9004) 

9004  FCPMAT< »i< ) 
WPITE(6,250)  PHIS 

250  FGPMAT('0','  PHI  STAR  =',2X,F10.7) 
WPITE(6,9005 )  N , RE Y, A • THO , VEL 

9005  FCPMATC  N  =',14,/,'  REY  =  '  ,  F  10  .  2  ,  6X,  •  AL  PHA  =  ', 

*  2F12.7,3X, '7HETA  = •  ,  F12 . 7T / t  ' VEL  =',F7.2) 
WRITE(6,9055)ASTAR 

9055  FGPMAT( '0'  ,' A-STAR  =',F12.7) 
C 

C  CALL  SUBROUTINE  TO  SOLVE  FOR  EIGENVALUES. 
C 

CALL  DE I  GEO { A , TH , R  E  Y , N , M 01 M , GR E , G I E , GRO , G 1 0  ,  XMAT , YMm T, 

*  8MAT,WVEC) 
C 

C  DETERMINE  WHICH  EIGENVALUE  IS  THE  lEAST  STABLE. 
C 

TEMP  =  -1D10 

MARK  =  1 
C 

CC  20  I=1,N 

IF(GRO(I )+AR*VEL.LT.TEMP)  GO  TO  20 

TEMP  =  GROU  )+AR*VEL 

ITEMP  =  I 
20  CCNTINUE 
C 

DO  40  I=ltN 

IFCGREC I )+AR*VEL.LT.TEMP)  GO  TO  40 

TEMP  =  GRE( I )+AR*VEL 

ITEMP  =  I 

MARK  =  2 
40  CCNTINUE 
C 

C  LIST  EIGENVALUES  FOR  ODD  E I G6NFUNCT I CNS 
C 

WRITE<6,9002) 
9002  FORMAT*///, 6X, 'GAMMA  REAL ', 5X,  ' GAMMA  IM AG '  ,  12X, • STAB ' ) 

WRITE(6,9006) 

9006  FORMAT* «0E  IGEN VALUES  FOR  ODD  EIGENVECTORS',/) 
C 

DC  50  1=1, N 
TEMP  =  GROU  )+AR*VEL 
WRITE* 6, 9000 )  GRO ( I  ) , GIO ( I ) ,TEMP 
I  FCI.EQ. ITEMP.AND.MARK.SQ.il  WRITE*6,90  01) 
5C  CCNTINUE 
90  30  FORMAT ( 'O' , 1P2D1 5. 4, 1PD20. 4 ) 
9001  FORMAT(  '  +  •  ,52X, •***»  ) 
C 

C  LIST  EIGENVALUES  FOR  EVEN  EIGENVECTORS. 
C 

WRITE (6, 9007) 

9007  FORMAT( • OEI GENVALUES  FOR  EVEN  EIGENVECTORS',/) 
DO  55  1=1, N 
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c 
c 
c 
c 


TEMP  =  GREU  )+AR*VEL 
*RITE(6*9000)  GRE(  I)»GIE(I)  tTEMP 
IMI..EQ.ITEMP.AND.MARK.EQ.2J  WRI  TE  (6,900  1 } 
55  CCNTINUE 


PUT  EIGENVALUES 
StePOUTINE  TO  DC? 


INTO  SINGLE  PRECISION  VECTORS 
PLOTTING  FOR  GDC  FUNCTIONS. 


TC    PASS    TO 


C 

c 
c 


CG    60    I=1,N 
GR4U  )    =    SNGL(GRO<  I  )  ) 
6C    GI4U )    =    SNGL(GIO(I  )  ) 
WRITE<6,9004) 
CALL    PLCTP(GR4fGI4fNt0) 
WRITE<6,9005)    Nt RE Y, A ,THD, VEL 
WRITE(6,9055  )     ASTAR 

SIMILARLY    PLOT    EIGENVALUES    FOR    EVEN     EI GENFUNCT I CNS 

DC    65    I=1,N 
GP4dJ     =    SNGL(GREd)) 
65    GI4d  )     =    SNGL(GIE(  I)  ) 
WRITE(6T9004) 
CALL    PL0TP(GR4,GI4,N,0) 
*RITE(o,9005  )    N, REY, AfTHD,VEL 
WRITE(6,9055)    ASTAR 

GC    TO    1 
130    WRITE(6,9004) 
STCP 
ENC 
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c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 


SUBROUTINE  DEIGEO 


PURPOSE 

DEIGEG  SOLVES  THE  LINEARIZED  N A  V IER-S TOKE S  EQUATION 
FOR  PCISEULLE  FLOW.   THE  INPUTS  TO  DEIGEO  ARE  THE 
STARRED  WAVE  NO.*  ALPHAt  THETA  *,  AND  STARRED 
REYNULDS  NUMBER. 
CEISEO  OUTPUTS  THE  EIGENVALUES  FOR  GAMMA. 

USAGE 

CALL  CEIGEO(ALPHA,THETA,REYNQ,N,MOIM,  rfREVENfWI EVEN, 
WROOD,WIOOO,CDM,DM,3M,WV) 

DESCRIPTION  CF  PARAMETERS 

THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM... 
ALPHA, THET At  REYNOTN, MDIM 

ALPHA  -  THE  PERTURBATION  WAVE  NUMBER  IN  THE  FLOW 
DIRECTION  (X).   (COMPLEXES) 

REYNO  -  THE  REYNOLDS  NUMBER  <REAL*8) 

N  -  THE  SIZE  OF  THE  MATRICES  WHICH  IS  EQUAL  T3 
(ND-D/2  WHERE  NO  IS  THE  NUMBER  OF  DIVISIONS 
ACROSS  THE  CHANNEL.   (NOTE...  DEIGEG  SOLVES  THE 
PR3BLEM  ACROSS  THE  HALF  CHANNEL  TWICE  -  ONCE  FCR 
ThE  EIGENVALUES  CORRESPONDING  TO  THE  EVEN 
EIGENFUNCTIONS  AND  ONCE  FOR  THOSE  CORRESPONDING 
TO  THE  ODD  EIGENFUNCTIONS. 

MDIM  -  THE  COLUMN  DIMENSION  OF  THE  MATRICES  WHICH 
DEI  GEO  USES.   MDIM  MUST  BE  .GE.  N 

THE  FOLLOWING  ARE  OUTPUT  BY  DEIGEO 
MREVENtWI EVEN tWROOOfW IOOO 

WREVEN,WIEVEN  -  THE  REAL  AND  IMAGINARY  PARTS  IF    "HE 
EIGENVALUES  CORRESPONDING  TO  THE  EVEN 
EIGENFUNCTIONS.   DIMENSIONED  TO  AT  LEAST  N. 
(REAL*8) 

WR0DD,WI0DD  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 
EIGENVALUES  CORRESPONDING  TO  THE  ODD 


EIGENFUNCTIONS. 
<REAL*8) 


DIMENSIONED  TO  AT  LEAST  N. 


THE  FOLLOWING  MATRICES  MUST  BE  INPUT  TO  DEIGEO  AS 
WORKSPACE. 

CCM(MOIM,MDI M)   (C0MPLEX*16) 
CM(MDIM,MDIM)   <REAL*8) 
BM<5,MDIM)   (C0MPLEX*16) 
WV(MDIM)   (CGMPLEX*lto) 

NCTES... 

THE  MATRICES  CAN  BE  OVERLAPPED  TO  CONSERVE  SPACE, 
FOR  EXAMPLE,  FOR  N  =  60... 

CCMPLEX*16  CDM(60,30,2) 
COMPLEX* 16  DM (60, 30), *V(6C) ,3M( 5,60) 
EQUIVALENCE  ( DM( 1 , 1 ) , COM ( 1 , 1 , 3 )  ), 
<BM<1, 1) ,C3M( 1,1,3) ) , 
WVEC(1),CDM(1,6,3> ) 

NOTE  THAT  IT  IS  ONLY  THE  ACTUAL  SIZE  OF  Tl-ESE 
WORKSPACES  THAT  IS  IMPORTANT,  NOT  THEIR  TYPE. 
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C  OTHER  ROUTINES  NEEDED 

C 

C  THE  FOLLOWING  ARE  CALLED  BY  CEIGED 

C 

C  CHM1E1,CHM2E2,MSET,CDMTIN,  BMSET, MULDBM , CSPLIT, 

C  EFESSCELRH1C 

c 

c 

c 

SUE ROUTINE  DE IGE 0 ( AL PHA, THETA, REYNO, N, MD I M , 
*     WREVEN,WIEVEN, WRODD,WlGCD, CDM,DM,3M, WV ) 

IMPLICIT  C0MPLEX*16(A-H,0-Z) 

DIMENSION  IVECQOOJ 

REAL* 8  WREVEN( 1),WIEVEN{1)TWRODD(1),WIOOD(1) 

REAL*8  CCM(1),DM( 1),8M<1) ,WV<1 i 

RE£L*3  REY,DELY, REYNO 

REAL*8  TH, THETA 

CCMMON  /  COEFNT  /  A  ,  TH , G , RE Y , DEL Y 

EXTERNAL  CHM1 El , CHM2E 1 
C 

C  THIS  SUBROUTINE  SOLVES  THE  EQUATION   YV  =  GXV   WHE^E 
C  X  AND  Y  ARE  MATRICES,  V  IS  THE  EIGENVECTOR  AND  G  IS  THE 
C  EIGENVALUE.   THE  EIGENVALUES  ARE  DETE~MI,\ED  ANC  PASSED 
C  BACK  TO  THE  CALLING  PROGRAM  IN  WRGDD,  wTDDD,  WREVEN  AND 
C  WIEVEN. 
C 

A  =  ALPHA 

TH  =  THETA 

REY  =  REYNO 
C 

C  SET  UP  MATRIX  X  FOR  ODO  EIGENVECTORS. 
C 

CALL  MSET(CDM,N,MDIM, 1,CHM2E1) 
C 

C  INVERT  MATRIX  X. 
C 

CALL  CDMTIN(N,CDM,MDI M,OETERM) 
C 

C  SET  UP  MATRIX  Y  IN  BAND  STORAGE  MODE  FOR  ODD  EIGENVECTORS. 
C 

CALL  3MSET(BM,N, MDIM,  1,CHM1E1) 
C 

C  MULTIFLY  MATRIX  Y  BY  THE  INVERSE  OF  MATRIX  X  TO  CONVERT 
C  TO  THE  STANDARD  EIGENVALUE  PROBLEM  WHICH  HAS  THE  FORM 
C  (Z-G)V  =  0   WHERE   Z  =  ( Y  )  ( INVE RSE ( X  )  )  . 
C 

CALL  MULDBMCCDM, 3M  ,N , 5  ,  MDI M, WV ) 
C 

C  SPLIT  MATRIX  INTO  REAL  AND  IMAGINARY  PARTS  ANC  CALL 
C  THE  SLEROUTINES  TO  FIND  THE  EIGENVALUES. 
C 

CALL  DSPLIT(N,MDIM,COM,CDM,DM) 

CALL  EHESSCCCDM, OMf lfNfNtMOIM, IVES) 

CALL  ELRH1C( CDM, DM,1  ,N ,N , MDI M, WR C CO, W IODC, I NE RR , I ER ) 

IF(INEPR.NE.O)  WRITE(6,9C00)  INERRTIER 
9O0C  FCRMAT( 'OERROR  NUMBER',17,'   ON  EIGENVALUE •,  I  7, /// ) 
C 

C  REPEAT  THE  SOLUTION  FOR  EIGENVALUES  FCR  THE  EVEN 
C  EIGENVECTORS 
C 

CALL  MSET(CDM,N,MDIM, 2,CHM2E1) 

CALL  CDMTIN(N,CDM,MDI M,DETERM) 

CALL  3MSET(8M,N, MDIM, 2,CHM1E1) 

C^LL    MULDBM(COM. BM,N, 5, MDIM, WV) 

CALL    DSPLIT(N,MDI M,CDM,CDM,DM ) 

CAL-L  -EHESSC(  COM,  DM,1  ,N,N,MDIM,  IVEC  J 

CALL  ELRH1C(CDM, DM , 1 , \, N , MDIM, WR EVEN , WI E VEN , I NERR , IER) 

IFUNERR.NE.O)  WRI  TE  (6  ,9000  )  INERR,IER 

RETURN 

ENC 
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C... SUBROUTINE    MULDBM 

C 

C 

C  PURPOSE 

C 

C        NULDBM  PERFORMS  THE  MATRIX  MULTIPLICATION  BETWEEN  A 

C        SQUARE  MATRIX  X  AND  A  BANDED  MATRIX  XB  WHICH  HAS 

C        BEEN  SET  UP  BY  SUBROUTINE  BMSET.   THE  RESULT  IS 

C        PLACEC  BACK  IN  X.   THAT  IS... 

C  X  =  (X)  (XB) 

c 

C      USAGE 

C 

C        CALL  MULDBM(X,XB,N,N3,MDIM,TEMPV) 

C 

C      DESCRIPTION  OF  PARAMETERS 

C 

C        THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM. 

C         NfMDIMtNBtXf XBfTEMPV 

C 

C        N  -  THE  SIZE  OF  THE  MATRICES. 

C 

C        MDIM  -  THE  DIMENSIONS  OF  THE  MATRICES  IN  THE  CALLING 

C  PROGRAM. 

C 

C        NB  -  THE  NUMBER  OF  BANDS  IN  THE  BANDED  MATRIX,  X3. 

C 

C        X  -  THE  SQUARE  N  BY  N  MATRIX.   CIMENSIGNED 

C  (MDIM, MDIM)  IN  THE  CALLING  PROGRAM   (COMPLEX* 16 ) 

C 

C        XB  -  THE  N  BY  N  MATRIX  WHICH  IS  STORED  IN  BANDEC 

C  FORM.   DIMENSIONED  (N3,MDIM)  IN  CALLING  PROGRAM. 

C  (C0VPLEX*16) 

C 

C        TEMP\/  -  A  VECTOR  WORKSPACE  WHICH  MUST  BE  PROVIDEC  BY 

C  THE  CALLING  PROGRAM,   DIMENSIONED  AT  LEAST  (N). 

C  (C0MP1_EX*16) 

C 

C        THE  FCLLCWING  IS  OUTPUT  BY  MULCBM. 

C 

C        X  -  SET  TO  THE  PRODUCT  OF  X  AND  XB  (MDIM,VDIM) 

C  <CGMPLEX*16) 

C 

C      REQUIRED  ROUTINES 

C 

C        NONE 

C 

c 

C.3 

c 

S L  BROUT I N E  MU LD3 M ( X , X  8 , N, W 3 , MD I M  ,  T  EM P  V ) 

COMPLEX* 16  X(MDIM,MDI M) ,X3(NB, MDIM) t TEMP V(MCIM), TEMP 

NehM  =  (NB-D/2 

NEHP  =  (NB+1J/2 
C 

C  LCCP  CVER  INDEX  I 
C 

CO  100  1=1, N 

c 

C  STCRE  COLUMN  I  OF  MATRIX  X  TEMPORARILY 
C 

CC  10  J=1,N 
10  TEMPV( J)  =  X(I,J) 

C 

C  FIND  PRODUCTS  FOR  FIRST  NBHM  SPECIAL  CASES,  THAT  IS  WHERE 

C  WHERE  THE  BANDED  MATRIX  DOES  NOT  HAVE  ITS  FULL  WIDTH 


C 


CO    22    J =1, NEHM 
TEMP    =     (ODOtODO) 
JJ    =    MBHM    +    J 
CC    21    K=1,JJ 
JJJ    =    JJ-K+1 
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c 
c 

c 

c 
c 


c 

c 

c 


21  TEMP  =  TEMP+TEMPV(K) *XB( JJJ,K) 

22  X(I,J)  =  TEMP 

COMPUTE  PRODUCTS  FOR  "REGULAR"  COMBINATIONS  CF  RO'rtS  ANO 
CGLUMNSt  ThAT  IS,  THOSE  THAT  ARE  NOT  TRUNCATED 
AT  THE  BEGINNING  OR  END  BY  THE  BOUNDARIES 

JF  =  N-NBHM 

CO  22  J=NBHPjJF 

TEMP  =  (OCCODO) 

CO  31  K=1,NB 

JJJ  =  NB-K+1 
21  TEMP  =  TEMP+TEMPV( J-NBHP+K ) *X3 ( J J J , J-NBHP+K ) 
32  X(  I, J)  =  TEMP 

FIND  PRODUCTS  FOR  LAST  NBHM  SPECIAL  CASES. 

DC  42  J=1,NBHM 
TEMP  =  (ODO.ODO) 
JJ  =  NB-J 
DO  41  K  =  1,JJ 

41  TEMP  =  TEMP+TEMPV(N-JJ+K)*XB(N3-K+1,N-JJ+KJ 

42  XUfN-NBHM+J)     =    TEMP 

100  CONTINUE 

RETURN 
ENC 
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C.  , SUBROUTINE  MSET 

C 

c 

C  PLPPOSE 

c 

C  MSET  SETS  UP  A  MATRIX  FOR  THE  FINITE  DIFFERENCE 

C  PROBLEM  OF  POISEULLE  FLOW  WITH  THE  PROPER  BOUNDARY 

C  CONDITIONS  FOR  THE  VELOCITY  VECTOR  POTENTIAL  FOR 

C  VISCOUS  FLOW. 

C 

C  USAGE 

C 

C  CALL  MSET(X,N, MDIM, MODE, CFMAT  ) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM 

C  N, MDIM, MODE. CFMAT 

C 

C  N  -  THE  SIZE  OF  THE  MATRIX.   IF  M3DE=0,  N  IS  EGUAL 

C  TO  THE  NUNBER  OF  POINTS  OF  THE  FINITE  DIFFERENCE 

C  MESH  ACROSS  THE  CHANNEL  NOT  INCLUDING  Tt-E  TWO 

C  BOUNDARIES.   IP  MGDE=1  OR  MODE=2,  N  IS  EQUAL  TO 

C  CNE-HALF  OF  THE  NUMBER  OF  INNER  POINTS  ACROSS  THE 

C  FULL  CHANNEL.   IN  THIS  CASE,  THE  CHANNEL  MUST  BE 

C  DIVIDED  INTO  AN  EVEN  NUMBER  CF  POINTS  SC  THAT  N 

C  WILL  BE    AN  INTEGER. 

C 

C  MDIM  -  THE  COLUMN  DIMENSION  QP  THE  MATRIX  X.   MDIM 

C  MUST  BE  .GE.  N. 

C 

C  MODE  -  IF  MODE=0,  THE  MATRIX  IS  SET  UD  FOR  THE  FULL 

C  CHANNEL.   IF  MODE=l  OR  MODE=2,  THE  MATRIX  IS  SET 

C  UP  FOR  THE  HALF  CHANNEL  AND  THE  BOUNDARY 

C  CONDITIONS  ARE  SET  SUCH  THAT  THE  ODD  OP  EVEN 

C  EIGENFUNCTIONS ,  RESPECTIVELY,  ARE  SOLVEC  FOR. 

C 

C  CFMAT  -  THE  NAME  OF  A  FUNCTION  SUBPROGRAM  WITH  TWO 

C  PARAMETERS,  K  AND  Y,  INDICATING  WHICH  TERM  JF  THE 

C  FINITE  DIFFERENCING  IS  DESIRED,  AND  THE  POSITION 

C  RELATIVE  TO  THE  CENTER  OF  the  CHANNEL.   MUST  BE 

C  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

C  (C0MPLEX*16) 

C 

C  THE  FCLLOWING  IS  OUTPUT  BY  MSET 

C 

C  X    -    THE    N    BY    N    MATRIX    INTO    WHICH    THE    COEFFICIENTS 

C  OF    THE    FINITE    DIFFERENCING    ARE    PUT.       MUST    BE 

C  DIMENSIONED     (MDIM, MDIM)     IN    THE    CABLING    PROGRAM. 

C  <C0MPLEX*i6) 

C 

C  NCTES... 

C 

C  THE  EIGENVALUES  AND  VECTORS  OBTAINED  3Y  USING  MSET 

C  TWICE,  WITH  MCDE=1  AND  MCDE=2,  ARE  THE 

C  SAME  AS  USING  MSET  ONCE  WITH  MCOE=0,  BUT  WITH 

C  N  TWICE  AS  LARGE. 

C 

C  OTHER  ROUTINES  NEEDED 

C 

C  NONE    -    EXCEPT    THE    FUNCTION    SUBPROGRAM    NAME    PASSED    IN 

C  THE    PARAMETER     'CFMAT'. 

C 

c 

c 

SLEROUTINE  MSET ( X, \,MD I M, MODE, CFMAT  ) 

REAL*8  REY,Y,OELY,DFLOAT 

CCNPLEX*lo  A,G 

REALMS  TH 

CCMMON  /  COEFNT  /  A,TH,G,REY ,DEL Y 

CCNPLEX*16   CFMAT 

COMPLEX*16  X(MDIM,MDIM) 
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c 

C  CCMPUTE  GRIC  SIZE  FOR  FINITE  DIFFERENCE  MESH  ACROSS 

C  HALF  CHANNEL  OR  FULL  CHANNEL 

C 

CELY  =  2D0/DFL3AT(N+1 ) 

IF(MC0E.EQ.1.0R.M0DE.EQ.2)  OELY  =  200/D FLC AT ( 2*N  + L ) 
C 

C  CHECK  IF  MATRIX  DIMENSIONED  LARGE  ENOUGH 
C 

IF(N.LE.MDIM)  GO  TO  1 

WPITE(6»9000) 
9000  FCPMATt'O*  *  *  ERROR  -  ARRAYS  NOT  DIMENSIONED  LARGE', 
*     •  ENOUGH  *  *  *' ) 

STOP 
C 

C  ZERO  ENTIRE  MATRIX 
C 

1  CC  10  1=1, N 

DC  10  J=1,N 
10  X(ItJ)  =  (ODO,ODO) 
C 

C  DC  SPECIAL  CASE  AT  DISTANCE  DELY  FROM  CHANNEL  WALL 
C  INCLUDING  6CUNDARY  CONDITIONS 
C 

Y  =  1D0-DELY 

X(l,l)     =    CFMAT(3,Y)+CFMATll ,Y) 

X(l,2)     =    CFMAT(4, Y) 

X(l,3)     =    CFMAT(5,Y) 
C 

C  DO  SPECIAL  CASE  AT  DISTANCE  2*DELY  FROM  CHANNEL  WALL 
C  INCLUCING  30UNCARY  CONDITIONS 
C 

Y  =  1D0-2D0*DELY 
X(2,l)  =  CFMAT(2,Y) 
X(2t2)  =  CFMAT(3,Y) 
X(2,3)  =  CFMAT<4, Y) 
X<2,4)  =  CFMAT(5,Y) 

C 

C  DO  ALL  REGULAR  POINTS  IN  BETWEEN,  THAT  IS,  THOSE  VALUES 

C  OF  Y  FCR  WHICH  ALL  5  FINITE  DIFFERENCE  GRID  POINTS  ARE 

C  INTERIOR  TO  THE  CHANNEL 

C 

IL  =  N-2 

DO  20  1=3, IL 

K  =  1-3 

Y  =    1Q0-DELY*DFL3AT(T ) 
DC    20    J=l,5 

20    X(I,K+J)     =    CFMAT ( J,Y) 
C 

C  FINALLY  DO  THE  TWO  SPECIAL  CASES  WHICH  OCCUR  EITHER  AT 
C  THE  CENTER  OF  THE  CHANNEL  DR.  AT  THE  OTHER  WALL,  DEPENDING 
C  ON  THE  VALUE  OF  MODE.   BOUNDARY  CONDITIONS  ARE  SET  UP 
C  DEFENDING  CN  MODE 
C 

Y  =  1D0-DELY*DFL0AT(N-1) 
X(N-l,N-3)  =  CFMAT(1,Y) 
X(N-i,N-2)  =  CFMAT(2,Y) 
X(N-lfN-l)  =  CFMAT(3,Y) 
X(N-1,N)  =  CFMAT(4,Y) 

IF(MODE.EQ.l)  X(N-1,N)  =  CFMAT ( 4  ,  Y ) -CFM AT ( 5 , Y ) 
IF(MQDE.EQ.2)  X(N-ltN)  =  CFMAT ( 4 , Y ) +CFM AT ( 5 , Y ) 
C 

Y  =  100-DELY*DFL0AT<N) 
X(N,N-2)  =  CFMAT( 1,Y) 
X(N,N-1 )  =  CFMAT(2,Y) 

IF(MrOE.EQ.l)  X(N,N-1)  =  CFMAT ( 2  ,  Y ) -CFM AT ( 5 , Y ) 

IF(M0DE.EQ.2)  X(N,N-1)  =  CFMAT ( 2  ,  Y ) +CFMAT ( 5  ,  Y ) 
X(N,N)  =  CFMAT(3,Y)+CFMAT(5,Y) 

IF(MODE.EQ.l)  X(N,N)  =  CFMAT ( 3 , Y ) -CFM AT < 4 , Y ) 

IF(M00E.EQ.2)  X(N,N)  =  CFMAT ( 3 , Y ) +CFMAT ( 4, Y ) 
C 

RETURN 
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ENC 
C SUBROUTINE    BMSET , 

C 

C 

C  PLPPOSE 

C 

C  THE    PURPOSE    OF    BMSET    IS    EXACTLY    THAT    OF    MSET    EXCEPT 

C  THAT    BMSET    TAKES    ADVANTAGE    OF    THE    BANDED    NATURE    OF 

C  THE    FINITE    DIFFERENCE    MATRICES    TO    CONSERVE    SPACE. 

C 

C  USAGE 

C 

C  CALL    BMSETU, N,MDIM, MODE, CFMAT  ) 

C 

C  DESCRIPTION    OF    PARAMETERS 

C 

C        THE  PARAMETERS  FOR  BMSET  ARE  THE  SAME  AS  THOSE  FOR 

C        MSET  WITH  THE  EXCEPTIuN  THAT  THE  MATRIX  X  MUST  BE 

C        DIMENSIONED  (5,MDIM)  IN  THE  CALLING  PROGRAM. 

C 

C      NCTE:   THE  PROCEDURE  IS  IDENTICAL  TO  THAT  OF  MSET. 

C        COMMENTS  HAVE  THEREFORE  NOT  BEEN  INCLUDED  IN  3MSET. 

C 

c 

c... 

c 

SUBROUTINE  BMSET ( X , N, MDIM, MODE , CFMAT  ) 
REAL*8  REY,Y,DELY,DFLCAT 
CCMPLEX*16   CFMAT 
CCVPLEX*16  AVG 
CCMPLEX*16  X(5,MDIM) 
REAL'S  TH 

CCMMGN  /  COEFNT  /  A ,TH , G, REY , DELY 
DELY  =  2D0/DFL0A7(N+1) 

IF(MODE.EQ.1.0R.MODE.EQ.2)  DELY  =  2D0/DFLCAT ( 2*N+ 1 ) 
IF(N.LE.MDIM)  GO  TO  1 
WRITE(6,900C) 
9000  FCRMATCO*  *  *  ERROR  -  ARRAYS  NOT  DIMENSIONED  LARGE', 
*     «  ENOUGH  *  *  *• ) 

STCP 
1  DC  10  1=1,5 
DC  10  J=1,MDIM 
10  X(ItJ)  =  (ODO,ODO) 

Y  =  1D0-DELY 

X(2,l)  =  CFMAT(3,Y)+CFMATQ  ,Y) 
X<4,1)  =  CFMAT(4, Y) 
X<5,1)  =  CFMAT(5, Y) 

Y  =  1D0-2D0*DELY 
X(2,2)  =  CFMAT(2,Y) 
X  (3,2)  =  CFMAT(3, Y) 
X(4,2)  =  CFMAT(4,Y) 
X(5,2)  =  CFMAT(5,Y) 
IL  =  N-2 

CO  20  1=3, IL 

Y  =    1D0-DELY*DFL0AT(I  ) 
DC    20    J =1,5 

20    X(J,I )  =    CFMAT( J, Y) 

Y  =    1D0-DELY*DFL0AT(N-1) 
X(  1,N-1 )  =    CFMAT(lt Y) 
X(2,N-1)  =    CFMAT(2, Y) 
X(3,N-1J  =    CFMAT(3,Y) 

X  (4,N-1)     =    CFMAT(4,Y) 

IF(M0DE.EQ.1  )     X(4,N-1)  =    CFMAT  (  4  ,  Y  ) -CFM  AT  (  5  ,  Y) 

IF(M0D£.EQ.2)     X(4,N-1)  =    CPMA T (4  ,  Y ) +CFM AT ( 5  ,  Y  ) 

Y  =    1D0-DELY*DFL3AT(N) 
X(1,N)  =    CFMAT(  1,Y) 
X(2,N)          =    CFMAT(2,Y) 

IF(MODE.EQ.l)  X(2,N)  =  CFMAT< 2 , Y ) -CFMAT ( 5 , Y ) 
IF(M0DE.EG.2)  X(2,N)  =  CFMAT( 2 , Y  ) +CFM AT < 5, Y  ) 
X(3,N)     =    CFMAT13, Y ) +C FMmT ( 5 , Y ) 

IF(MODE.EQ.l)  X(3,N>  =  CFMAT ( 3, Y  ) -CFMAT ( 4, Y  ) 
IF(M0DE.EQ.2)     X(3,N)     =    CFMAT ( 3 , Y  ) +CFMAT ( 4 , Y  ) 
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RETURN 
ENC 
C • SUBROUTINE    DSP  LIT 

C 
C 
C      PURPOSE 

c 

C        DSPLIT  TAKES  A  MATRIX  OF  CCMPLEX*16  NUMBERS  AND 

C        SPLITS  IT  INTO  TWO  MATRICES,  ONE  CONTAINING  THE  REAL 

C        PART  CF  THE  ORIGINAL  MATRIX,  ANC  ONE  CONTAINING  THE 

C        IMAGINARY  PART. 

C 

C      USAGE 

C 

C        CALL  CSPLIT(N,MDIM,A,AREAL, AIMAG) 

C 

C      DESCRIPTION  OF  PARAMETERS 

C 

C        N  -  THE  SIZE  OF  THE  MATRIX  A,  AN  N  BY  N  SCUARE 

C  MATRIX. 

C 

C        MDIM  -  THE  COLUMN  DIMENSION  OF  MATRIX  A 

C 

C        A  -  THE  INPUT  MATRIX.   MUST  BE  DIMENSIONED  MDIM  BY 

C  AT  LEAST  N  IN  THE  CALLING  PRCGRAM   <CGMFLEX*16) 

C 

C  AREAl,AIMAG    -    THE    OUTPUT    MATRICES    CONTAINING    THE 

C  REAL  AND  IMAGINARY  PARTS,  RESPECTIVELY,  OF 

C  MATRIX  A.   MUST  3E  DIMENSIONED  (MDIM, MDIM)  IN  THE 

C  CALLING  PROGRAM. 

C 

C      NCTES... 

C 

C        MATRIX  A  AND  MATRIX  AREAL  MAY  OVERLAP  IF  THEY  ARE 

C  DIMENSIONED  IN  THE  CALLING  PROGRAM  AS  FCLLGWS... 

C 

C  C0MPLEX*16  A(MDIM,MDIM) 

C  REAL*8    AREAL (MCIM, MDIM)  ,AIMAG( MDIM, MCIM) 

C  EQUIVALENCE( A(l,l) ,AREAL(1, 1) ) 

C 

C      OTHER  ROUTINES  NEEDED 

C 

C        NONE 

C 

C 

C... 

r 

SUBROUTINE  CSPLI T ( N , MDIM , A , AR , AI  ) 

REAL* 8  A<2»MDIM, MDIM), ARC MDIM TMOI*)fAI (MDIM, MDIM) 

C 

CO  1  J=1,N 

DC  1  1=1, N 

AR(I,J)  =  A(l, I, J  ) 
1  Aid, J)  =  A(2,  I,  J) 
C 

RETURN 

ENC 
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c, 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c, 

c 


c 
c 

c 


c 
c 

c 


SLBROUTINE  CDMTIN  (CATEGORY  F-l J 
PLRPOSE 

INVERT  A  C0MPLEX*16  MATRIX 
LSAGE 

CALL  CDMTIN<N»AfNDIM,DETERMl 

DESCRIPTION  OF  PARAMETERS 

N       -  ORDER  OF  C0MPLEX*16  MATRIX  TO  3E  INVERTED 
(INTEGER)  MAXIMUM  »N'  IS  100 

A       -  C0MPLEX*16  INPUT  MATRIX  (DESTROYED).   THE 
INVERSE  OF  «A'  IS  RETURNED  IN  ITS  PLACE 

NDIM    -  THE  SIZE  TO  WHICH  «A'  IS  DIMENSIONED 

(ROW  DIMENSION  OF  'A*  ACTUALLY  APPEARING 
IN  THE  DIMENSION  STATMENT  OF  USER'S 
CALLING  PROGRAM) 

DETERM  -  C0MPLEX*16  VALUE  OF  DETERMINANT  OF  'A' 
RETURNED  BY  CDMTIN. 

REMARKS 

MATRIX  'A1  MUST  8E  A  CGMPLEX*8  GENERAL  MATRIX 

IF  MATRIX  'A*  IS  SINGULAR  THAT  MESSAGE  IS  PRINTEO 

•N«  MUST  BE  .LE.  NDIM 

SLBRCUTINES  ANO  FUNCTIONS  REQUIREC 

ONLY  BUILT-IN  FORTRAN  FUNCTIONS 

METHOD 

GAUSSIAN  ELIMINATION  WITH  COLUMN  PIVOTING  IS  USED. 
THE  DETERMINANT  IS  ALSO  CALCULATED.  A  DETERMINANT 

OF  ZERO  INDICATES  THAT  MATRIX  'A«  IS 

SINGULAR. 


SLeROUTINE  CDMTIN  ( N, A ,NDI M , DETERM ) 
IMPLICIT  REAL*8  (A-H.O-Z) 

COMPLEX* 16  A (NDIM, NO  I M ) , P I VOT ( 100  )  , AMAX , T , S *A P , 
*     DETERM. U 
INTEGER *4  I  PI VOT ( 100) ,INDEX( 100 r 2  ) 
REAL*3  TEMP,ALPHA(100) 

INITIALIZATION 

CETERM  =  (lDOtODO) 
DC  20  J  =  1,N 
ALPHA(J)  =  ODO 


CC  10  1=1. N 
10  ALFHA( J)=ALPHA( J)+A(J ,I)*DCCNJG( A(J, I) ) 

20 


ALPHA ( J)=DSQRT( ALPHA ( J) ) 
IFIVOT( J)=C 
CO  600  1=1, N 


SEARCH  FOR  PIVOT  ELEMENT 

AMAX  =  (OCO,ODO) 
DO  105  J=1,N 

IF  (IPIVOT(J)-l)  60,105,60 
60  CO  100  K=1,N 

IF  (  IPIVOT(K)-l)  80,100,740 
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80    T=yP=AMAX*DCQNJG(AMAX  )-A( J , K ) *DCGNJG( A( J  ,  K)  ) 

IF(TEMP)  85,  85,  100 
8  5  IB0W»J 
ICCLUM=K 
AMAX=A( J,K) 
100  CCNTINUE 
105  CONTINUE 

IPIV0T(IC0LUM)=IPIV0T(IC0LUM)+1 
C 

C       INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 
C 

IF(IROW-ICCLUM)  140,  260,  140 
140  DETERM=DETERM 
CO  200  L=1,N 
SWAP=A(  IROWtU 
A(IROW,L)=A< ICOLUM,L) 
200  A(  ICOLUM,L)=SWAP 
SWAP=ALPHA(IRQW) 
ALPHA(IROW)=ALPHA(ICOLUM) 
ALPHA(iCOLUM)«SWAP 
260  INCEX(  I,1)  =  IR0W 

INCEX(I,2)=IC0LUM 
PIVOT(Il*A(ICOLUM,ICOLUMI 
L=PIVOT(I ) 

TEMP=PI VOT(I )*0C3NJG(PIV0T( I ) ) 
IF(TEMP)  330,  720,  330 
C 

C       CIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 
C 

330  A(ICGLUM,ICOLUM)  =  (1C0,0D0) 
DC  350  L=1,N 
L=PIVOT(I ) 
350  A(ICOLUM,L)=MICOLUM,L)/U 
C 

C       REDUCE  NDN-PIV3T  ROWS 
C 

38C  CO  550  L1=1,N 

IF(Ll-ICOLUM)  400,  550,  400 
400  T=A(L1, ICOLUM) 

A(L1,ICCLUM)  =  (OCO,ODO) 
CO  450  L=1,N 
U=A(ICOLUM,L) 
450  A(L1,L)=A(L1,L)-J*T 
550  CCNTINUE 
600  CONTINUE 
C 

C       INTERCHANGE  COLUMNS 
C 

620  DC  710  1=1, N 
L=N+1-I 

IF(INDEX(L,1)-INDEX(L,2) )  630,  71C,  630 
630  JR0W=INDEX(L,1) 

JC0LUM=INDEX(L,2 ) 
CO  705  K=1,N 
SWAP=A(K, JRCW) 
A(K, JROW)=A(K, JCOLUM) 
A(K,JCOLUM)=SWAP 
705  CONTINUE 
710  CCNTINUE 

RETURN 
720  WRITE(6,730) 

730  F0RMAT(20H   MATRIX  IS  SINGULAR) 
74C  RETURN 
END 
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C.EHESSC. 0 , 

c 

C            FUNCTION  REDUCTION    OF    A    COMPLEX    MATRIX    TO 

C  UPPER    HESSENBERG    FORM, 

C          USAGE  -    CALL    EHESSC ( AR , AI . K , L , N  .  I  A, I D ) 

C          PARAMETERS         AR  -    INPUT/OUTPUT    MATRIX    OF     CIMENSION 

C  N    BY    N. 

C  ON  INPUT  CONTAINS  THE  REAL 

C  COMPONENTS  OF  THE  REDUCED  HESSEN 

C  BERG  FORM  IN  UPPER  TRIANGULAR 

C  PORTION  ANC  THE  DETAILS  OF  THE 

C  REDUCTION  IN  LOWER  TRIANGULAR 

C  PORTION. 

C                  AI  INPUT/OUTPUT  N  BY  N  MATRIX 

C  CONTAINING  THE  IMAGINARY  COUNTER 

C  PARTS  TO  AR,  ABOVE. 

C                  K  -  INPUT  SCALAR  CONTAINING  THE  ROW 

C  AND  COLUMN  INDEX  OF  THE  STARTING 

C  ELEMENT  TO  BE  REDUCED  BY  ROW 

C  SCALING.  FOR  UNBALANCED 

C  MATRICES  SET   L  =  N. 

C                 L  INPUT  SCALAR  CONTAINING  THE  ROW 

C  AND  COLUMN  INDEX  OF  THE  LAST 

C  ELEMENT  TO  BE  REDUCED  BY  ROW 

C  SCALING.  FOR  UNBALANCED 

C  MATRICES  SET  K  =  1. 

C                 N  INPUT  SCALAR  CONTAINING  THE  ORDER 

C  OF  THE  MATRIX  TO  BE  REDUCED. 

C                  IA  INPUT  SCALAR  CONTAINING  ROW 

C  DIMENSION  CF  AR  AND  AI  IN  THE 

C  CALLING  PROGRAM. 

C                 ID  OUTPUT  VECTOR  OF  LENGTH  L  CONTAIN 

C  ING  DETAILS  OF  THE 

C  TRANSFORMATION. 

C    PRECISION  -  SINGLE/DOUBLE 

C    CODE  RESPONSIBILITY  -  T.J.  AIRD/E.W.  CHOJ 

C    LANGUAGE  -  FORTRAN 

C 

C    LATEST  REVISION  -  FEBRUARY  7,  1S73 
C 

SUBROUTINE  EHESSC  ( AR , AI , K, L , N , I  A , ID ) 
C 

DIMENSION  AP(IA,l)tAI(IA,l)»ID(l),71(2),T2(2) 

DOUBLE  PRECISION    AR , AI , XR , XI , YR , YI , T 1 , T 2 , ZERO 

CCMPLFX:Stl6  X  Y 

ECU I  VALENCE  (X,T1<1),XR),(T1(2)»XI),(Y,T2(1),YR), 

1  (T2(2),YI) 

CATA  ZERO/O. OCO/ 

LA=L-1 

KP1=K+1 

IF    (LA    .LT.    KP1)     GO    TO    45 

DO    40    M=KPltLA 
I  =  M 

XR=ZERO 
XI=ZERO 
DO    5    J=MrL 

IF(DA83UR(  J,M-1  ))+DABS(AI<  J  T  M- 1 )  )  .LE. 


* 

DA6S(XR)+DABS(XI  ) 

1 

GO    TO    5 

XR=AR( Jf M-l ) 

XI=AI( J,M-1 ) 
I=J 
CONTINUE 

ID(M)=I 

IF    (  I     .EQ.    M)     GO    TO    20 

INTERCHANGE    ROWS    ANC    COLUMNS 

ARRAYS    AR    ANC    AI 

MM1=M-1 

DO    10    J=MM1,N 

YR=AR( IT J  1 

ARdtJ  )=AR(M,J  ) 

AR(M,J)=YR 
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YI=AI(I, J) 
AI  (I  ,J)=AI  (M,  J) 
AKM,  J)  =  YI 
10     CONTINUE 

DO  15  J  =  1,L 
YP=AR( J, I ) 
AR(  J, I  )=AR( J»M) 
AR( J,W)=YR 
YI=AI( J, I) 
AI( J, I  )  =  AI ( J,M) 
AI( J,M)=YI 
15     CONTINUE 

END    INTERCHANGE 
20  IF    (XR    .EQ.    ZERO    .AND.    XI     .EQ.    ZERO)     GO    TO    40 

MP1=M+1 
DO    35    I=MP1,L 
YR=AR(I,M-1 ) 
YI=AI(I,M-1) 

IF    (YR    .EQ.     ZERO    .AND.    YI     .EQ.    ZERO)     GO    TO    35 
Y  =  Y/X 

AR(I  ,M-1)=YR 
AI(I,M-1)=YI 
CC    25    J=M,N 

AR{  I, J)=AR(I , J)-YR*AR (M, J )+YI* A  I ( M , J ) 
AH  It J)=AI( I, J)-YR*AI(M,J)-YI*AR(M,J) 
25        CONTINUE 

DC  30  J=1,L 

AR{ J,M)=AR(Jf M)+YR*AR(J, I ) -YI *  A  I < J  ,  I  ) 
AI(J,M)=AI(J,M)+YR*AI(JTI)+YI*AR(JTI) 
20        CONTINUE 
25     CONTINUE 
40  CONTINUE 
45  RETURN 
END 
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C.ELPH1C D. 

C 

C    FUNCTION  -  COMPUTATION  OF  ALL  EIGENVALUES  GF 

C  A  COMPLEX  UPPER  HESSENBERG 

C  MATRIX. 

C    USAGE  -  CALL  ELRH1C ( HR ,HI , K , L , N , I H, WR , Wl , 

C  INFER, IER) 

C    PARAMETER     HR  -  INPUT  MATRIX  OF  DIMENSION  N  BY  N 

C  CONTAINING  THE  REAL  COMPONENTS 

C  OF  THE  COMPLEX  HESSENBERG  MATRIX 

C  HR  IS  DESTROYEQ  ON  CUTPJT. 

C                  HI  INPUT  MATRIX  OF  DIMENSION  N  BY  N 

C  CONTAINING  THE  IMAGINARY  COUNTER 

C  PARTS  TO  HE,  ABOVE.   hi  IS 

C  DESTROYED  ON  OUTPUT. 

C                  K  INPUT  SCALAR  CONTAINING  THE  LOWER 

C  BOUNDARY  INDEX  FOR  THE  INPUT 

C  MATRIX.   FOR  UNBALANCED  MATRICES 

C  SET  K  =  1. 

C                  K  INPUT  SCALAR  CONTAINING  THE  UPPER 

C  BOUNDARY  INCEX  FOR  THE  INPUT 

C  MATRIX.   FOR  UNBALANCED  MATRICES 

C  SET  L  =  N. 

C                 N  INPUT  SCALAR  CONTAINING  THE  ORCER 

C  OF  THE  MATRIX. 

C                  IH  INPUT  SCALAR  CONTAINING  THE  ROW 

C  DIMENSION  OF  MATRICES  HR  AND  Hi, 

C  IN  THE  CALLING  PROGRAM. 

C                 WR  OUTPUT  VECTOR  OF  LENGTH  N  C3NTAIN 

C  ING  COMPONENTS  OF  THE 

C  EIGENVALUES. 

C                  WI  OUTPUT  VECTOR  OF  LENGTH  N  CONTAIN 

C  ING  IMAGINARY  COMPONENTS  OF  THE 

C  EIGENVALUES. 

C                  INFER  -  OUTPUT  SCALAR  CONTAINING  THE  INDEX 

C  OF  THE  EIGENVALUE  WHICH 

C  GENERATED  THE  TERMINAL  ERROR. 

C  N*l  INDICATED  THE  EIGENVALUE 

C  RECORDED  IN  THE  OUTPUT 

C  PARAMETER,  INFER,  CDULD 

C  NOT  3E  DETERMINED  AFTER  30 

C  ITERATIONS.  IF  THE  J-TH 

C  EIGENVALUE  COULD  NOT  BE  SO 

C  DETERMINED,  THEN  THE  EIGEN 

C  VALUES  J  +  l, J+2,  ...,N 

C  SHOULD    BE    CORRECT. 

C          PRECISION  -    SINGLE/DOUBLE 

C         REQ'D    IMSL    ROUTINES    -    UERTST 

C         CODE    RESPONSIBILITY    -    T.J.    AIRD/E.w.    CHOU 

C          LANGUAGE  -    FORTRAN 

C 

C    LATEST  REVISION  -  MARCH  22,     1973 
C 

SUBROUTINE  ELRH1C  ( HR , HI  ,  K, L , N , I H , wR , Wl  ,  INF ER , I E R ) 
C 

DIMENSION  HR(IH,1),HI(IH,1),WR{1),WI(1) 

DIMENSION  TK2)  ,T2(  2)  ,T3(2) 

CCMPLEX*16  X,Y,Z 

DOUBLE  PRECISION  HR , HI , WR , W I ,EPS , ZR, ZI , T i , T 2, T3 , RDEL P 

DOUBLE  PRECISIGN  ZERO ,ONE ,T WO , SR ,S I , XR, X  I , YR , RI , TR , TI 

EQUIVALENCE  ( X, Tl ( 1 ) ,XR ) ,  (Tl ( 2 ) , X  I  ) , 

1  (Y,T2( 1)  ,Y«)  ,  (T2(2) ,YI) , 

2  (Z,T3(1)  ,ZR) ,  (T3(2),ZI ) 
DATA  TWO/2. 000/ 

CAT  A    ZERO,CNE,RDELP/O.DO, 1. DO , Z3410000 000000000/ 

INFER=0 

IER  =  0 

DO  5  1=1, N 

IF  (I  .GE.  K  .AND.  I  .LE.  «_ )  GO  TO  5 

WR(I  )=HR(I, I ) 

Hid  )=HI(I,I) 
5  CONTINUE 

79 


10 


15 


NN 
TR 
TI 

IF 
IT 

N^ 
IF 


NF 
DC 


=  L 

=  ZERG 
=  ZERQ 


(NN 

s=o 

1=NN-1 
(NN  .EQ.  K) 


SEARCH  FCR  NEXT  EIGENVALUE 
LT.  K)  G3  TO  9005 


20 
25 
30 


35 
40 


45 


C 
C 


50 
55 

60 


CO 
M  = 
IF 
IF 

IF 
SR 
SI 
XR 
XI 
IF 
YR 
YI 
Z= 
IF 
X  = 
SR 
SI 
GC 
SR 
SI 
DC 


CC 
TR 
TI 
IT 


L-NN+K 

20  LL=K,NM1 

M=NPL-LL 

MM1=M-1 

IF  (DABS(HR(M 
RDELP*(DA8S( 
DA8S(HR(M, 
NTINUE 
K 

(M  .EQ.  NN)  G 

(ITS  .EQ.  30) 


.0  TO  25 

LOOK  FOR  SINGLE  SMALL  SUB-DIAGONAL 
ELEMENTS 


,MM1M  +  DABS  (HI  (M,MM1)  )  .LE. 
HR(MMlfMMD)  ♦  DABSiHl  (MM1,MM1)  ) 
M))  +  DABS(HI(M,M) )) )  Gu  TO  30 


(ITS  .EQ. 
=HR(NN,.\IN) 
■HUNNvNN) 

=HR(NM1,NN)*HR 
=HR(NM1,NN)*HI 

(XR  .EQ.  ZERO 
=  (HR(NM1,NMD- 
=(HI(NM1,NM1)- 
CDSQRT(DCMPLX( 

( YR*ZR+YI*ZI 
X/(Y+Z) 
=SR-XR 
=SI-XI 

TO    40 
=DABS(HR(NN,NM 
=DABS(HI(NN,NM 

45    I*KtNN 

HR(I  ,1  )=HR(  I, 

Hid, I  )=HI(I, 
NTINUE 
=TR+SR 
=TI+SI 
S=ITS+1 


0  TO  110 
GO  TO  115 

FORM  SHIFT 
10  .OR.  ITS  .EQ.  20)  GO  TO  35 


(NN,NM1)-HI (NH1,NN)*HI(NN,NV1 ) 
(NN,NMl)+HI(NMl,NNJ*HR(iMN,Nyi) 

.AND.    XI    .EQ.     ZERO)    GO    TO    40 
SR)/TWO 
SI  I /TWO 

YR**2-YI**2+XR, rWO*YR*YI*XI ) ) 
.LT.     ZERO)     Z=-Z 


1  )  )+DABS(HR(NMl,NN-2)  ) 
1  )  )+DABS(Hl  (NMl,i\N-2)  ) 


XR 
YR 

ZR 
NN 
IF 

DO 


=DABS(HR(NM1,N 
=DABS (HR(NN,NM 
=DABS <HR(NN,NN 
J=NM1-M 
(NMJ     .EQ.    0) 


I  )-SR 
D-SI 


LOOK    FOR    TWO    CONSECUTIVE     SMALL 

SUB-DIAG1NAL    ELEMENTS 
Ml  )  )+DABS(HKNMl,NMl  )  ) 
1  J  )+DABS(Hl  (NN,NM1I ) 
)  )+DABS(HI(NN,NN)  ) 

GO    TO  55 

FOR  MM=NN-1  STEP  -1  UNTIL  M+l  DO 


CC 
Mf 


50  J=1,NMJ 

MM=NN-J 

M1=MM-1 

YI=*YR 

YR=DA8S(HR(MM,M1)  ) +D ABS ( HI ( MM , Ml)  ) 

XI=ZR 

ZR-  XR 

XR=DABS(HR(M1 

IF  (YR.LE.RDE 
NTINUE 
=  M 


,M1  )  )+DABS(HKMl,Ml)  ) 
LP*ZR/YI*(ZR+XR+XI) )  GO  TO  60 


MP1=MM+1 

DG  85  I=MP1,NN 
IM1=I-1 
XR=HR(  IM1, IM1 
XI=HI( IM1,IM1 
YR=HR(If IM1) 
YI=HI( I,IM1) 


TRIANGULAR  OECCMPOS IT  ION 


80 


IF(DA3S(XR)+DABS(XI ).GE.CABS( YR)+DA3S(YI  )  )  GO  TO  70 
C  INTERCHANGE  ROWS  OF  HR     AND  HI 

DO  65  J=IM1,UN 
ZR=HR(IM1, J) 
HR(IM1, J)=HR(I, J) 
HR(  I,J)=ZR 
ZI=HI(IM1, J) 
HI( IM1, J)=HI (I, J ) 
HI(I,J)=ZI 
65     CONTINUE 
Z=X/Y 
WR(I )=ONE 
GO  TC  75 
7C     Z=Y/X 

WR(I  )=-ONE 
75     HR(I,IM1)=ZR 
Hl(ItIMl)*ZI 
DO  80  J=I,NN 

HR( I ,J)=HR( I ,J) -ZR*HR(IM1, J)+ZI*HI (IM1, J) 
HI(I,J)=HI ( I, J)-ZR*HI (IMit J)-ZI*HRUM1,  J) 
80     CONTINUE 
85  CONTINUE 
C  COMPOSITION 

DO  105  J=MPltNN 
JM1=J-1 
XR=HR( J,JM1) 
XI=HI( J, JM1) 
HR(JTJM1)=ZER0 
HI( J,JM1)=ZER0 
C  INTERCHANGE  COLUMNS  OF  HR  AND  HI  IF 

C  NECESSARY 

IF  (WR(J)  .LE.  ZERO)  GO  TO  95 
DO  90  I=MrJ 

ZR  =  HR(I, JM1  ) 
HR( I ,JM1 )=HR(If J ) 
HR(I,J)=ZR 
ZI  =  HI(  If  JMU 
HI(I,JM1)=HI (I, J  J 
HI(I,J)=ZI 
90     CONTINUE 
95     DO  100  I=M,J 

HR(I,JM1)  =  HR(  I,  JM1)+XR*HR(  I  ,  J)-XI  *HI  (  I ,  J  )■ 
HI(  I  ,JM1)  =  HI  (I,  JM1)+XP.*HI  (  I  ,J)+XI*HP(  I, J) 
100     CONTINUE 
105  CONTINUE 
GC  TO  15 
C  A  RCCT  FOUND 

110  WR(NN)=HR(NN,NN)+TR 
WI(NN)=HI(NN,NN)+TI 
NN=NM1 
GC  TO  10 
C  SET  ERROR-NO  CONVERGENCE  TC  AN 

C  EIGENVALUE  AFTER  30  ITERATIONS 

115  INFER=NN 
IER=129 
900C  CONTINUE 

CALL  UERTST  ( I ER , 6HELRH1C ) 
9005  RETURN 
END 
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C  SUBROUTINE  UERTST  (IER,NAME) 

C 

C-UEPTST LIBRARY 

C 

C  FUNCTION             -  ERROR  MESSAGE  GENERATION 

C  USAGE                -  CALL  UERTST (  I  ER , NAME ) 

C  PARAMETERS    I ER     -  ERROR  PARAMETER.  TYPE  +  N   WHERE 

C  TYPE=  128  IMPLIES  TERMINAL  ERROR 

C  64  IMPLIES  WARNING  WITH  FIX 

C  32  IMPLIES  WARNING 

C  N      ERROR  CODE  RELEVANT  TO 

C  CALLING  ROUTINE. 

C  NAME    -  INPUT  VECTOR  CONTAINING  THE  NAME 

C  OF  THE  CALLING  ROUTINE  AS  A  SIX 

C  CHARACTER  LITERAL  STRING. 

C  LANGUAGE             -  FORTRAN 

C 

C  LATEST  REVISION      -  JANUARY  18,  1974 


C 


SLEROUTINE  UERTST ( I ER, NAME ) 

DIMENSION  ITYP(5T4) ,IBIT(4) 

INTEGER*2  NAMEOl 

INTEGER  WARN, WARF, TERM, PRINTR 

ECU  I  VALENCE  ( I  BIT ( 1 ) , WARN ) ,<IBIT(2),WARF),(IBIT(3), 

*  TERM) 

CATA       ITYP       /'WARNS'ING    •,'  '»*  ,f<  »i 

*  'WARN' ,' ING( •,  'WITH'  , 'FIX'  .'  )     ', 

*  'TERM' , • INAL* ♦ ■      ■ , '     ' , •     ■ , 

*  IBIT       /  32,64,128,01 
DATA      PRINTR     /  6/ 


IER2=IER 

IF  (IER2  .GE.  WARN)  GO  TO  5 

IER1=4 
GC  TO  20 
5   IF  (IER2  .LT.  TERM)  GO  TO  10 

IER1=3 
GO  TO  20 
10   IF  (IER2  .LT.  WARF)  GC  TO  15 


NON-DEFINED 


ERMINAL 


WARNING(WITH  FIX) 


IER1=2 
GC  TO  20 
C  WARNING 

15   IER1=1 
C  EXTRACT  »N» 

20   IER2=IER2-I3IT(IER1) 
C  PRINT  ERROR  MESSAGE 

WRITE  (PRINTR, 25)  ( IT YP ( I , I ER 1 ) , I =1 , 5 ) , N AME  ,  I ER2, I ER 
25   FORMAT(»  ***  I  M  S  L(UERTST)  ***   •  ,5 A4 , 4X , 2A2, 4X , I  2 , 
*    •  (IER  =  ',  13,  ■  )• ) 
RETURN 
END 
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C FUNCTION    CHM1E1    AND    CHM2E1 

C***#" (CART  ESI  AN*  COORDINATES^ 

C 
C 

C  PURPOSE 

C 

C  CHM1E1  AND  CHM2E1  RETURN  THE  VALUES  (C0MPLEX*16)  OF 

C  THE  COEFFICIENTS  FOR  THE  MATRICES  IN  THE  FINITE 

C  DIFFERENCE  FORM  OF  THE  LINEARIZEO  NAV I ER-STOKES 

C  EQUATION  FOR  POISEUILLE  FLOW.   BOTH  FUNCTIONS  RESULT 

C  FROM  THE  LINEAR  COMBINATION  OF  EQUATION  1  AND 

C  EQUATION  3  TO  ELIMINATE  THE  VELOCITY  VECTCR 

C  POTENTIAL  COMPONENT  G  AND  ARBITRARILY  SETTING  THE 

C  COMPONENT  F  TO  ZERO.   Su ,  THEY  APE  THE  COEFFICIENTS 

C  FOR  THE  VECTOR  POTENTIAL  COMPONENT  H.   CHM2E1 

C  RETURNS  THE  TERMS  WHICH  ARE  COEFFICIENTS  CF  THE 

C  EIGENVALUE,  GAMMA,  AND  CHM1E1  RETURNS  THE  REMAINING 

C  TERMS, 

C 

C  USAGE 

C 

C  XI  =  CHM1E1(K,Y) 

C  X2  =  CHM2E1(K,Y) 

C 

C  (CHM1E1  AND  CHM2E1  MUST  BE  DECLARED  C0MPLEX*16  IN 

C  CALLING  PROGRAM) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  PARAMETERS  MUST  BE  SET  BY  THE  CALLING 

C  PROGRAM 

C 

C  K  -  INDICATES  THE  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH  RELATIVE  TO  THE  CENTRAL  PGINT  IN  THE  CENTRAL 

C  DIFFERENCING  SCHEME.   IF  THE  DIFFERENCE  IS  BEING 

C  FORMED  ABOUT  THE  N-TH  POINT  THEN  K=l  REFERS  TO 

C  THE  POINT  N-2,  K=2  REFERS  TO  THE  POINT  N-l, 

C  K  =  3  REFERS  TO  N,  K=4  REFERS  Tu  N  +  l,  AND  K=5 

C  REFERS  TO  N+2. 

C  K  -  INDICATES  WHICH  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH  IS  REFERRED  TO  THE  CENTRAL  POINT.   IF  THE 

C  DIFFERENCE  IS  BEING  FORMED  ABOUT  THE  N-TH  POINT 

C  THEN  K=l  REFERS  TO  THE  POINT  N-2,  K=2  REFERS  TO 

C  THE  POINT  N-l,  K=3  REFERS  TO  N,  K=4  REFERS  TC  N+l, 

C  AND  K=5  REFERS  TO  N+2. 

C 

C  Y  -  THE  VALUE  OF  THE  POSITION  RELATIVE  TO  THE  CENTER 

C  OF  THE  CHANNEL.   THE  TWO  BOUNDARIES  ORE  AT  Y=+l 

C  ANC  Y=-l. 

C 

C  OTHER  ROUTINES  NEEDED 

C 

C  NONE 

C 

C 

c 

c 


FUNCTION  CHM1EKK.Y) 

IMPLICIT  C0MPLEX*16  (A-H,0-Z) 

CCMMCN  /  COEFNT  /  A, TH, G, REY ,OEL 

REAL*8  REY, Y, DEL 

REAL*8  TH,DUR 
C 

C  THE  FGLLOWING  FUNCTIONS  (Ml)  EVALUATE  THE  COEFFICIENTS 
C  OF  THE  DERIVATIVES  OF  H  FOR  ALL  TERMS  EXCEPT  THCSE 
C  CONTAINING  GAMMA. 


C 


CH4MKY)     =  A*EI/REY 

CH2MKY)     =  -1.5DO*A**2*EI2*(ldO-Y**2)+2DO*AEI*U**2)/ 
*                 REY 

CHOMKY)    =  -AEI*<  (A**2)*(1.5DO*AEI*(lDO-Y**2)-(A**2) 
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*  /REY)+3D0*AEI ) 
C 

C  THE  REGAINING  FUNCTIONS  (M2)  EVALUATE  THE  COEFFICIENTS 
C  OF  THE  DERIVATIVES  OF  H  WHICH  ARE  ALSO  COEFFICIENTS 
C  OF  THE  EIGENVALUE  GAMMA. 
C 

CF2M2(Y)  =  AEI 

CFCM2(Y)  =  AEI  *  A**2 

CUR  =  0.0 

DL  =  DCMPLX(DUR,TH) 

EI  =  CDEXP(DU) 

AEI  =  A*EI 

EI2  =  C0EXP(2*DU) 
C 

C  SET  UF  THE  FINITE  DIFFERENCE  VALUES  FOR  INDEX  K  FOR  Ml 
C 

GC  TO  ( 1,2,3,2,1),K 

1  CHM1E1  =  CH4M1(Y)/DEL**4 
GC  TO  100 

2  CHM1E1  =  -4D0*CH4M1(Y)/DEL**4-H:H2M1(Y)/DEL**2 
GC  TO  100 

2  CHM1E1  =  6D0*CH4Mim/DEL**4-2DC*CH2Ml(Y)/DEL**2 

*  +CHOMKY) 
100  RETURN 

C 

C  SET  UP  THE  FINITE  DIFFERENCE  VALUES  FOR  INDEX  K  FOR  M2 

C 

ENTRY  CHM2E1(K,Y) 

CLR  =  0.0 

CL  =  OCMPLX(DURtTH) 

EI  =  CDEXP(DU) 

AEI  =  A*EI 

EI2  =  CCEXP(2.0*DU) 

GC  TO  (11,12,13, 12,11),K 

11  CHM2E1  =  (ODO,ODO) 
GC  TO  200 

12  CF^2E1  =  CH2M2(Y)/DEL**2 
GC  TO  200 

12  CHM2E1  =  -2C0*CH2M2(Y)/DEL**2+CHOM2(Y) 
200  RETURN 
EMC 
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